VERTICAL CURRENT STRUCTURE IN THE GREAT LAKES Vincent E. Noble Joseph C. Huang James H. Say lor Final Report Federal Water Pollution Control Administration Grant WP-OIO67 Special Report No, 57 of the Great Lakes Research Division The University of Michigan Ann Arbor^ Michigan 1968 TABLE OF CONTENTS Page LIST OF TABLES v LIST OF FIGURES vi INTRODUCTION 1 I. CHARACTERIZATION OF THE CIRCULATION DYNAMICS OF LAKE MICHIGAN (Vincent E. Noble) k II. CURRENT METER STUDIES OFF SHEBOYGAN, WISCONSIN, AUGUST 196? (Vincent E. Noble) 12 III. DETERMINATION OF EDDY VISCOSITY AND EDDY DIFFUSIVITY IN LAKE MICHIGAN (Joseph C. Huang) 26 Introduction 26 The Magnitude of Eddy Viscosity and Eddy Diffusivity in Lake Michigan 28 An estimate of eddy viscosity by Reynolds stress and vertical velocity gradients 28 An estimate of the upper bound value of eddy viscosity by the result from statistic theory 29 An estimation of vertical eddy viscosity from current measurement near the bottom 29 An approximation of the vertical eddy viscosity by surface wind 51 A measurement of the eddy diffusivity by the dispersion of marked particles 3^ Instantaneous sources 55 Continuous sources 55 A measurement of the horizontal diffusion by drogue studies 58 Comparison of the Values of Eddy Viscosity and Eddy Diffusiv- ity with Other Studies 58 Conclusion ij-5 IV. NUMERICAL EVALUATION OF STERN'S CIRCULATION MODEL (James. H. Say- lor) k6 Development of Stern *s Model k6 Application of Stern's Model to Lake Michigan 5I Effect of density difference, Ap 55 Effect of T 56 Effect of depth of thermocline 56 iii TABLE OF CONTENTS (Concluded) Page V. THERMAL CURRENT STRUCTURE IN LAKE MICHIGAI^, ASSOCIATED WITH THE SPRING WARM-UP SEASON— A THEORETICAL STUDY (Joseph C. Huang) 59 Introduction 59 Modeling Considerations and Assumptions 6o Mathematical Formulation 6k Nondimensional equations 65 Solutions of the Zeroth p- Order Approximation 68 Temperature 68 Boundary layer scaling 68 Interior solutions 69 Boundary-layer solutions 71 Results of the Zeroth p-Order Solution 7^ Temperature 7^ Velocity 85 The First p- Order Solutions 86 Conclusions 91 Future Plan 92 IV LIST OF TABLES Tables Page 5-1 « Eddy viscosity on two current-meter buoy stations (17-2^^- August 1967 at 8 miles off Sheboygan). 50 5-2 « Estimated vertical eddy viscosity by current measurement near bottom, 51 5-5* Estimated vertical viscosity by wind speed, 52 5-4„ General information and effective eddy diffusivity of dye patch studies. 5U 5-5« General information and effective eddy diffusivity of continu- ous source dye plume study. 59 5-6. General information and eddy diffusivity of drogues, study in Lake Michigan. k2 h"!^ Calculated eddy diameters, phase velocities, and internal wave periods of upper layer in Lake Michigan. 55 U-2. Calculated eddy diameters, phase velocities, and internal wave periods for upper layer in Lake Michigan with density difference changed. 55 ^-5< Calculated eddy diameters, phase velocities, and internal wave periods for upper layer with wind stress changed. 56 k-k. Calculated eddy diameters, phase velocities, and internal wave periods in Lake Michigan with the depth of thermocline changed. 56 5-1. Nomenclatures and nondimensional parameters. 65 5-2. Local coordinates in different regions. 70 V LIST OF FIGURES Figure Page 1-1. Surface temperature contours from Airborne Radiation Thermometer flight of 18 October I966. 5 1-2. Surface temperatures and 0-10 m dynamic height currents from BT transects taken by R/v INLAND SEAS, October I966. 6 1-5. Surface temperature contours from Airborne Radiation Thermometer flight of 25 October 1967. 8 1-U. Surface temperatures and 0-10 m dynamic height currents from BT transects taken by R/V INLAND SEAS^ 26 October I967. 9 1-5, Temperature structure of a vertical cross-section of Lake Michi- gan corresponding to Airborne Radiation Thermometer flight tracks. 10 2-la Location of current meter stations, showing contours of bottom topography. I5 2-2. Comparison of Milwaukee (MKE) winds and i|-l/2 hour current vec- tor averages for stations I and II at 0-10 m depths (15 August - 2k August 1967). 15 2-5. Dynamic height current components and corresponding ^-l/j hour current meter averages, I7 August I968. I7 2-U. Drogue observations and quarter-period current meter averages, 18 August 1968. 18 2-5. Drogue observations (showing drogue displacement and net trans- port velocity), dynamic height current components, and corre- sponding ^-1/5 hour current meter averages. I9 2-6. Dynamic height current components and corresponding ^-l/5 hour current meter averages^ 21 August I967. 21 2-7. Dynamic height current components and ^-l/5 hour current meter averages, 23 August I967. 22 5-1. A typical distribution of dye concentration across a dye patch. ^k VI LIST OF FIGURES (Concluded) Figure Page 5-2. Lateral growth of the dye patch with time. 55 5-5* Longitudinal growth of the dye patch with time. 5^ 5-^. Continuous source apparatus. 57 5-5* Lateral growth of dye plume at a constant distance from the source. 40 5-6. Lateral growth of dye plume with distance from source. kl i+-l. Schematic two-layer model. 51 U-2. Specific parameters of the two-layer model. 5^ 5-1. East-west cross- section of the central portion of Lake Michigan. 6l 5-2. Cross- section of a long, symmetrical trapezoidal Lake Michigan model. 62 5-5- Boundary regions of the Lake Michigan model. 70 5-^. Surface temperature patterns during spring warming of Lake Michigan from available data, 195^-68. 75 5-5^ Vertical temperature structure of central portion of Lake Michigan during warming period, I965-68. 79 5-6. Generalized horizontal zonal currents. 85 5-7» Meridional and vertical circulation. 86 5^8. Surface currents in Lake Michigan, 29 June 1955 (after Ayers et al. 1958). 87 5-9» Brogue studies for surface current measurements near eastern shore . 88 5-10. Drogue studies for surface current measurements near western shore . 90 5-11. Circulation associated with "Thermal Bar" (after Rodgers I965). 91 VI 1 INTRODUCTION The research program carried out under this research grant (WP-OIO67) was a two-year continuation of the studies initiated under a previous U.S. Public Health Service grant (WP-0079^). The final report describing the previous work was prepared by Noble (I966). The previous report shows the failure of Ekman- type models in synthesizing current meter records from buoy stations when using the wind measured at the station as an input stress. The previous study further suggests that geostrophic forces defined by the temperature-dependent density structure of Lake Michigan may be the dominant factor in determining the three- dimensional circulation patterns of the lake. Four major field experiments were carried out (in cooperation with other projects and individuals) under this grant. The two most significant operations were thermal mapping flights carried out by the Ant i- Submarine Warfare Environ- mental Prediction Services (ASWEPS) of the U.S. Naval Oceanographic Office. The results of the l8-21 October I966 flight program of the ASWEPS Super Constella- tion^ EL COYOTE (showing the existence of a major circulation feature^ the Coyote Current) have been described in detail in a report by Nobel^ £t al. (1967). The results of the second flighty on 25 October 1967^ are compared with those of the 1966 flight in a paper in preparation by V. E. Noble and J. C. Wilkerson. A third major field experiment was carried out in August 1967j> in coopera- tion with Dr. G. E. Birchfield of Northwestern University. A summary of the results of this study was presented by Birchfield and Noble. -^ In this study^ two current meter stations were established in Lake Michigan for a two-week period^ and comparisons were made between the current meter records, drogue measurements, and dynamic height current calculations. Spectral analyses and inertial period averages were also computed from the current meter records. The fourth field experiment consisted of a series of measurements of the temperature, eddy viscosity, eddy diffusivity, and current structures in the nearshore edge of the lake during the spring warming period of I968. These measurements were carried out by J. C. Huang. The results will be used as part of the basic input data and parameters for his theoretical investigations of lake circulation, and will be compared with the results obtained in connection with another project by Noble and Anderson (1968) during the spring of I967. In addition to the experimental program, two theoretical studies have been carried out in connection with this grant by students of the Department of Statistical analysis of currents at two nearby stations in Lake Michigan, summer 1967.. G. E. Birchfield and V. E. Noble. Paper presented at 11th Conf . on Great Lakes Res., Milwaukee, I968. Meteorology and Oceanography of The University of Michigan. The first was the output from a Special Problems' course in which James H. Saylor calculated numerical results, based on Lake Michigan data, from Stern's (1965) model for the interaction of a uniform wind stress with a geostrophic vortex. The second is a Ph.D. dissertation being written by Joseph C. Huang. Mr. Huang is doing a theoretical analysis of a simplified Lake Michigan model using the early spring heat input to the lake as the prime current- gene rating force. Early re- sults from this model (Huang's dissertation is not expected to be completed until the summer of I969) indicate that the geostrophic currents resulting from the spring temperature structure of the lake may be as great as previous compu- tations for wind-driven currents. Those results correspond with observations, and with the concepts within Stern's (I965) model. In the proposal for this grant, it was anticipated that the data from a special triangular array of current meter stations set in Lake Huron during the winter of I965-66 would provide a basis for investigation of the Eulerian- Lagrangian transformation inherent in the extrapolation of current meter data, and of correlations between data from pairs of current meters. Instrumentation problems caused a lack of data from this source, and therefore it was necessary to rely upon measurements made under this grant for subsequent analysis and in- terpretation. References NOBLE, VINCENT E. I966. VERTICAL CURRENT STRUCTURE IN TEE GREAT LAKES • Spec. Rep. No. 27, Great Lakes Res. Div., University of Michigan, k2 p. NOBLE, VINCENT E., et al. I967. NAVOCEANO FLIGHT OVER LAKE MICHIGAN, 17-21 October I966. Spec. Rep. No. 5I, Great Lakes Res. Div., University of Michigan, 17 p. + 2k figs. NOBLE, VINCENT E. and R. F. ANDERSON. I968. Temperature and current .structure in the Grand Haven, Michigan, vicinity during thermal bar conditions. Proc. 11th Conf. on Great Lakes Res., Internat. Assoc, for Great Lakes Res. In press. STERN, MELVIN E. I965. Interaction of a uniform wind stress with a geostrophic vortex. Deep-Sea Research, 12:555-567. I. CHAEACTERIZATION OF THE CIRCULATION DYNAMICS OF LAKE MICHIGAN Vincent E. Noble The results of previous investigations carried out in connection with this program tend to indicate that the circulation dynamics of Lake Michigan and, by inference, the other Great Lakes, are dominated by geostrophic forces rather than the wind stresses during the summer season. Analysis of temperature records from the Great Lakes-Illinois River Basins Project thermographs installed at buoy stations in Lake Michigan during the winter season of 1965-6^ indicated that the prevailing westerly wind drift caused an eastward surface current on Lake Michigan, a sinking on the east side of the lake basin, a westward bottom current, and an upwelling on the west side of the basin (Heap and Noble I966) . In the spring, the "thermal bar" (Rodgers I965; Noble and Anderson I968), which is a consequence of the spring heating of the lake, defines a general circula- tion around the edge of the lake basin. During the onset of the summer season, with the development of the thermocline, the wind stress-geostrophic stress interactions become extremely complex. Critical testing of theoretical models depends heavily upon our ability to make adequate density measurements in both time and space. Analysis of data from single current meter stations has failed to demonstrate a direct rela- tionship between wind and current. Not enough is yet known about the spatial structure of currents to interpret the relationship between data from current meter stations separated by distances of 5 to 10 miles. Surface temperature patterns obtained by synoptic observations from aircraft imply certain types of circulation features when interpreted with respect to dynamic height current computations made from BT transects. However, this interpretation is subject to errors inasmuch as the dynamic height computations depend on the essential assumption that the current field is defined by the density (temperature) field, and therefore it is only self-consistent that the dynamic height currents should show a strong correlation with the surface temperature structure. To the extent that dynamic height computations show a good qualitative correlation with cur- rents measured with drogues and current meters, the fundamental assumption of the dominance of the geostrophic forces is valid. The question yet to be re- solved is the interaction between the wind and geostrophic forces. The persis- tence of small details in the surface temperature structure of the lake for periods of three weeks (Heap and Noble I966; Noble I967) give credence to the postulate of dominance of the lake circulation dynamics by geostrophic forces. Figure 1-1 shows isotherms contoured from 5-mile averages of surface temperature along 25 east-west flight lines spaced at 10-mile intervals from the south end of Lake Michigan on 18 October I966. Subsequent flight operations showed that the general temperature pattern persisted for at least three subse- quent days. Figure 1-2 shows the dynamic height currents computed for the 0-10 m 42*00'-|- 8e»oo' Figure 1-1. Sairface temperature contours from Airborne Radiation Thermometer flight of 18 October I966. 20 10 10 I 87» 35' &?• 15' I 86« 55' I 86» 35* INLAND SEAS 18 OCTOBER 1966 43* 57.2' 87* 35' 87»I5 86» 55' INLAND SEAS 20 OCTOBER 1966 43* 57.2* — 44* 00* 87* 35' 87* 15' 86» 35* INLAND SEAS 19 OCTOBER 1966 Figure 1-2. Surface temperatures and 0-10 m dynamic height currents from BT transects taken by R/v INLAND SEAS, October I966. layer, compared with surface temperature readings taken from BT transects along the 45 '^ 57, 2* N flight line on 18, 19, and 20 October I966 by R/v INIAp SEAS. These dynamic height currents indicated a strong northward current parallel to the strong temperature gradient on the west side of the northern half of the lake. This circulation feature was also indicated by a distinct visual water mass separation along the gradient. The water was warm and green to the east of the gradient, and cold and blue to the west. The circulation feature ^as named the Coyote Current in honor of the NAVOCEANO research aircraft, the Super Constellation EL COYOTE. If the circulation of Lake Michigan is dominated by geostrophic (temperature- dependent density) forces, it would be expeated that the Coyote Current would appear each fall in the northern end of Lake Michigan. Therefore another flight was scheduled during the fall of I967 to see if the current did recur. The early fall period of I967 was much colder than that of 1966, and the temperature mapping flight could not be scheduled until 25 October. Therefore the 1967 flight was equivalent to being two weeks later in the season than that of 1966. The surface temperature contours from the I967 Coyote flight are shown in Figure I-5. Due to heavy weather, the INLAKD SEAS could not make a BT transeci on the day of the flight. Figure l-k shows 0-10 m dynamic height currents and the surface temperature profile from the INLAM) SEAS BT transect at ky^^J.^^'N on 26 October I967. The 1967 temperature gradient was stronger (of the order of ^''C). The temperature gradient turned offshore near Sheboygan, Wis., in I966, and more nearly paralleled the shoreline in I967. Figure 1-5 shows the temperature struc- ture of a vertical section of the lake as determined by the BT transects of 1966 and 1967- Because the development of the surface temperature structure was not followed continuously through the fall and early winter periods of 1966 and 1967^ it is not possible to make definitive judgments based on the compari- son of measurements made on 18 October I966 and 25 October I967. However, the appearance of the strong temperature gradient feature covering a large portion of the north- south extent of the lake basin, accompanied by a water mass separa- tion as evidenced by color changes visible from the aircraft, gives a strong indication that a similar phenomenon was observed in both years. Additional small-scale experiments were planned and executed to provide experimental foundations for testing and evaluation of theoretical models that would lead to an understanding of the circulation dynamics of the lake. These experiments were continued observations of the current regime associated with the early spring temperature structure, and a 2-week period of measurement using two current meter stations and comparison current drogue and BT observations. Figure 1-5 • S-urface temperature contoiirs from Airborne Radiation Thermometer flight of 25 October I967. 10 1 11 'I '■ ' 10 43«'57.3' INLAND SEAS 26 OCTOBER 1967 Figure 1-k. Surface temperatures and 0-10 m dynamic height currents from BT transects taken by r/v INLAND SEAS, 26 October 196?. (a) INLAND SEAS 18 OCTOBER 1966 go 70 QO 90 IJO -50 5^<C^ -100 -150 ==;^ \ ^^ -7=^ ';'C^^ 2 ^'^'^^^u***^'' 87«30' I 8 7" 00' I (b) INLAND SEAS 20 OCTOBER 1966 go go 90 ||0 86<'30 I 43''57.3 (c) INLAND SEAS 26 OCTOBER 1967 503070 e^S* , 10" 9" 9° 87*'30' 87*00' 86<»30 Figure 1-5. Temperature structure of a vertical cross-section of Lake Michigan corresponding to Airborne Radiation Thermometer flight tracks, (a) 18 October 1966. (b) 20 October 1966. (c) 26 October I967. 10 References HEAP, J. A. and V- E. NOBLE. I966. GROWTH OF ICE ON LAKE MICHIGAN. Spec. Rep. No. 26, Great Lakes Res. Div., University of Michigan, 9^ p. NOBLE, V. E. 1967- Evidences of geostrophically defined circulation in Lake Michigan. Proc. 10th Conf. Great Lakes Res., Internat. Assoc, for Great Lakes Res., p. 289-298. NOBLE, V. E. and R. F. ANDERSON. I968. Temperature and current structure in the Grand Haven, Michigan, vicinity during thermal bar conditions. Proc. 11th Conf. on Great Lakes Res., Internat. Assoc, for Great Lakes Res. Jn press. RODGERS, G. K. I965. The thermal bar in the Laurentian Great Lakes. Univer- sity of Michigan, Great Lakes Res. Div. Pub. No. 15:558-563. 11 II. CURRENT METER STUDIES OFF SHEBOYGAN, WISCONSIN, AUGUST I967 Vincent E, Noble During the period of I6 August through 2k August IS&J , two current meter buoy stations were established in 280 ft of water about 8 miles offshore from Sheboygan, Wis. The station locations and relative bottom topography in that area of the lake are shown in Figure 2-1. The buoy stations were set a half mile apart in roughly an east-west direction with current meters set nominally at the surface and at 10-m depths. The actual current speed rotors were at 2 and 12 m, respectively. Both current meter stations were obtained on loan from the Federal Water Pollution Control Administration, Great Lakes- Illinois River Basins Project. The experiment was designed to compare current drogue obser- vations and dynamic height current computations with current meter readings, and to provide an insight into the distance through which current meter readings can be extrapolated and into the Eulerian-Lagrangian transformation represented by current drogue and current meter measurements. Station I was set at l^i-i^O EST on 16 August I967 at a position of kykh.VE and 87''32.0'W. The station depth was 2&\- ft. Station II was set at I530 on 16 August 1967 at a position hyhk.2''^ and 87^52. 6'W. The water depth at sta- tion II was 280 ft. The station position was one-half mile at bearing 287*^ from station I. The movie film on which the data from the surface meter at station I were recorded was prematurely exposed, therefore, no useful data were obtained for the first ^1-5 hours after the station was installed. The other three current meters operated satisfactorily throughout the full period. The buoy station configuration was similar to that used by the Great Lakes- Illinois River Basins Project. In this instance, however, the subsurface buoy (which is used to provide tension on the meter string) was set so that it was just awash at the surface of the water. The thermocline was rather shallow in the region of measurement so that the 12-m current rotor was just about at the thermocline, generally being in the metalimnion. The weather during the period of measurement was generally unfavorable so that far less comparison data were obtained than had been desired. On I7 August only BT's were taken for dynamic height current computations. On the l8th only drogue observations were made. On the 19th the ship was weathered-in at the dock. On the 20th both drogue and BT observations were made. On the 21st, BT*s only. The ship was laid in on the 22d, with only BT*s again on the 25d. The stations were recovered on the 2^th. The winds were generally southwest to west on the l6th and 17th, with a frontal passage on the l8th producing northwest winds through the 19th. The wind went back to southwest to west on the 20th with a second cold front passage with north to northeast winds on the 21st remaining northeast through the 22d and going southeast on the 25d. 12 -4^50 87'50' SHEBOYGAN ■43M0 Figure 2-1. Location of current meter stations, showing contours of bottom topography. 15 The full statistical analyses of the current meter data have been discussed by Birchfield and Noble, The statistical analyses^ performed by Dr. G. E* Birchfield, Department of Engineering Sciences, Northwestern University, include determination of the mean transport current and a measure of the dispersion of the instantaneous current vectors during full and quarter- inert ial periods (17-I/5 hours, and ^-l/5 hours, respectively). In this report, the quarter-period means are compared with drogue measurements, dynamic height current computations, and wind records from Milwaukee and Muskegon. The wind records from the land stations are used to document the general wind regime and significant frontal passages. Figure 2-2 compares the wind records from the land stations with the quarter- period averages from each of the current meters. The quarter-period current meter averages discussed in this section of the report are the result of a first analysis carried out by Dr. Birchfield. The quarter-period averages were computed by taking the times that the current meters were set as t = 0. Therefore, the 4-l/5-hour averages for stations I and II are displaced in time by 50 minutes (the elapsed time between the setting of stations I and II) . The current meter averages are being calculated for coin- cident periods, and a complete report of the statistical analysis of the current meter data and drogue comparisons is being prepared by Birchfield and Noble. The complete report will discuss the current meter averages and vector disper- sions for the coincident periods, and will contain spectral analyses of the current meter data. The winds at Milwaukee (MKE) were S from I5 August through 0000 I8 August. A frontal passage occurred around O5OO on I8 August, with winds shifting to N at 15 kts (ships on the lake reported 18-25 kts, N) . The MKE winds diminished on the 19th and went light, SW on the 20th (a ship reported JO i^"ts at 190° at 2100 on 20 August). MKE winds were W-SW, 5-10 kts through 1200 21 August, with a strong frontal passage around I5OO on 21 August (ships reported 15-55 kts at 060°-020° behind the front during the period 0900-2100 on 21 August). MKE winds went east, decreasing through I80O on 2k August. (Ship reports gave 21 kts at 030° for 2100 on 22 August, and 09 kts at 010° for 2100 on 25 August.) The ship reports used for this study were carferry reports taken approximately 50 miles NE of the buoy stations. As shown in Figure 2-2, the quasi- steady currents between the fronts do not show a simple relation to the wind direction. The response of the current meter records to a frontal passage is the development of a rotational characteristic lasting for a period of I8 hours. (The rotational characteristic develops after a time lag of the order of 12 hours.) It is suggested that the complex relation . between the current vectors and wind field may be explained in Ekman transport and acceleration of geostrophic eddies as described by Stern ^s (I965) model. 2 See footnote 1, p. 1. Ik Figure 2-2 Comparison of Milwaukee (MKE) winds and 4-1/2 hour current vector averages for stations I and II at 0-10 m depths (15 August - 24 August 1967) Figures 2-5 through 2-7 show the results of the comparative measurements between the drogues, dynamic height current computations, and the quarter iner- tial period averages of the cu!i:rent meter records. The quarter period averages give the average current vector determined from a 4-l/5-hour interval. The dynamic height current computations yield currents between pairs of BT*s taken within a time interval of the order of l/2 hour within the comparison current meter period. The net drogue movements indicated by the dotted lines on the figures are used to determine a net transport vector. Figure 2-5 shows the relation between dynamic height currents and quarter period averages of the current meter records for 17 August 1967* For station I the averaging interval was O80O to 1220 hours and for station II, 085O to 1510 on 17 August (period 5) The surface current data for station I were lost due to premature exposure of the film but records are available for 10 m at station I, and and 10 m at station II. The average currents from stations I and II are 2O-5I cm/sec. Dynamic height currents are computed from three BT casts taken between 0850 and 0952 on I7 August. The dynamic height computations give only the component of current speeds normal to the path between BT casts. For this comparison period the components determined by dynamic heights were k.9 and 5*6 cm/sec, respectively, for the average current in the to 10-m level between the three BT casts. The dynamic height currents were calculated by the method of Ayers (1957), and assume a reference depth of 70 m. The computed speed components were within 20^ of those shown by the average values for the current meter stations. The dynamic height component computed for a position 2 miles east of the stations showed good agreement with the current meter average, but the component computed for a position I-I/2 miles west of the stations showed a direction I80** from that given by the station data. The ship winds at 0815 on I7 August were I9 kts, with direction 195 °« Figure 2-k shows a comparison between drogue observations and quarter period current meter averages for 18 August 1967. The four current drogues were set between 0955 and IO55 and recovered between 1^^-15 and 15^1-1 on I8 August. A frontal passage had occurred Just before setting the drogues, and only two positions of the drogues were determined — the set point and the recovery point. Drogue speeds were generally to the north with a net transport velocity ranging from 5-7 to 21.5 cm/sec. For station I a quarter period average encompassed the time be- tween 1000 and 1^20 on I8 August with surface current of 19*5 cm/sec and 10-m current of 55*6 cm/sec (period 11). The surface current speed was 5O.8 cm/sec with 10-m speed of 52.8 cm/sec at station II. The current drogues consisted of ^-ft square sheet metal current crosses suspended from modified net stake floats. The center of the current cross was at about the same depth as the current speed rotor for the surface meters. The drogue showed a good correlation with both speed and direction with the surface current meter of station II but very poor correlation with the direction of station I. The winds were generally north and at lUoO the ship wind was north 16 ^ ^ GO 00 a> 0) c> ^- o 0) E K 3 O Z 1 1 UJ UJ (o cc: o: o H LlI UJ < X K h- -J o UJ UJ CL S :e CO UJ X o o z u o LiJ -J Z z UJ UJ o UJ UJ 3 3 :s q: q: O o < cr q: o o z , 3 3 q: cr >" H O O o o O GO ^ u o CVJ to 00 bD •\ a w •H ■p rd ^ a o c^ ft H CQ TiJ ^1 fl U •H o :S o 'd ^ a 00 ai On CQ H -P a -P (D GQ o gj f < o o H •^ a ^ (D CQ ^ cu ^ bD , o ^ o P a) -p \> CQ ^ 5 "^^v..^ bD S •H Sh o (D ^ -P a ^ •H O ^ 4^ -d a SU (U 0) ft S ^ ^ CQ -p o a • 0) CVJ g ^ o (U CVl •\ S H o ir\ •H 1 CA F^ -=t H rO 00 IT LEGEND O METER - Surfoce METER - 10 meter CURRENT CURRENT DROGUE DROGUE DISPLACEMENT DYNAMIC HEIGHTS BT 32.8 43* 46' 19.3 43044, 1/2 87» 34' miles 43042. 87" 30* Figure 2-k. Drogue observations and quarter- period current meter averages, 18 August I968. The net drogue displacements and average transport velocities are shovra for a 5-hour interval. Wind N, 19- 25 kts. 18 ^ o I I (T a: ill UJ I- H UJ UJ UJ o < -J a. CO <0 o ro o fO a z Ui o 111 -J Z Z UJ UJ UJ UJ 3 ID a: q: o o (T cr o o 3 3 q: (T o o Q Q 2 < z O CQ CM ro ^ r t 00 CM • e GO O -fO 00 O H I 00 CM fO o 00 (M H •H is: O o •H (D -P +5 CQ o CM O H I O ON O 19 CD o o ro CVJ o ro 0) 0) o i^- o o k- £ H 3 O 2 (f) \ 1 lU :^ UJ C/) QC oc o h- UJ UJ < X h- H -J o UJ UJ a. UJ X Q Q LlI LU _J 2 H Z LU LU o liJ LU 3 3 s q: q: O CD < cr q: o O z 3 3 q: q: >■ K O O Q Q Q CD ilro CvJ CO c\J ■ o -P CM CO TIJ a •H ^ 6 • ro (U 1^ CD o •H • ^d '5' (U j^ 0) nd PM CQ H O 4J CJ d O <u O ^ CVJ ^ ITv ro o 1 o CM CX) ON (D H P«4 -P CQ O OJ o ^ H ro On o H 00 S I r t 20 LEGEND CURRENT METER - Surface CURRENT METER - 10 meter T5- DROGUE , DROGUE DISPLACEMENT > DYNAMIC HEIGHTS O BT o 43''46' 2.7 18.3 3.9 0\.^ |2.4 _ o 87°34' 1/2 — miles O 87° 32' 43° 42' 87'>30' Figure 2-6. Dynamic height current components and corresponding i^-l/5 hour current meter averages, 21 August I967. Current speed in cm/sec. Wind S, 7 kts. 21 LEGEND CURRENT METER - Surface CURRENT METER - 10 meter -T' DROGUE — X DROGUE DISPLACEMENT a* DYNAMIC HEIGHTS O BT 51.4 / 25.9 21.3 4.2 \= 43°46' O 43044. 33.7 1/2 • miles 86° 34' 43° 42' 86°30' Figure 2-7. Dynamic height cairrent components and U-I/5 hour current meter averages, 23 August I967. The dynamic height coraponents were computed for a time at the end of averaging period 57 (solid vectors) and the beginning of period 58 (broken vectors). Wind SE, 15 kts. 22 19-25 kts. It is to be noted that all four drogues, and the surface current averaged at station II for the quarter period, were all to the north, that is, into the wind. Figures 2-5^ and b show a comparison between drogues, dynamic height cur- rents, and current meter averages for 20 August 1967* The five drogues were set on an east-west line through the current meter stations between 0912 and 0935* The drogues were fixed between 1525 and l4^0 and recovered between I655 sind 17^0. Figure 2- 5a shows the first half of the measurements between the set and the first fix of the drogues, and Figure 2-5b shows the second half between the fix and recovery. Figure 2- 5a gives the quarter period average currents for the time period between 09^0 and 1^00 for station I (period 22) and between 1050 and 1^50 for station II (period 22). The dynamic height current component for a position 1 mile west of the current meter station shows good agreement (speed: 0.7 cm/sec) for the northerly component with the net transport of the drogue set at the same location (8.6 cm/sec to the north of east). The speed component (to the north at 0.1 cm/sec) computed for a position I-I/2 miles to the east of the stations is in the proper sense, but does not agree as well with its corresponding net drogue transport (9.8 cm/sec to the northeast). Of particular interest is the net south transport of the drogue set 2 miles to the west of the stations. There is good agreement between the drogues set at the current meter sites and the ^-l/j-hour meter averages. Figure 2- 5b compares the quarter period current meter records from 1^00 to 1820 for station I (period 25) and 1^-50 to I9IO for station II (period 25) with the net transport currents given by the motion of the drogue between their fix and recovery points. These two direct current measurements are compared with the to 10-m current calculated by the dynamic height method. In this case, there is less correlation between the surface current meter averages, and be- tween these averages and the drogue transports and the dynamic height current components. (Only two dynamic height components are calculated for this inter- val, shown as 2.5 and 0.8 cm/sec to the east.) Figure 2-6 gives the results of the experiments on 21 August I967 comparing the 10-m dynamic height currents with the quarter period average from station I in the interval of 0720 to 11^0 (period 27) and the interval of 08IO to I25O at station II (period 27). Figure 2-7 gives the results for 25 August in which dynamic height currents for the to 10-m layer are compared with current meter averages for two adjacent U-hour periods. The BT casts for the computation of the dynamic height currents were taken at the end of period 57 in the averaging sequence, therefore, the two adjacent periods 57 and 58 are compared with the dynamic height currents. For station I the two intervals are 02J^0 to 07OO and 0700 to 1120, and for station II O55O to 0750 and O75O to 1210 on 25 August. Of interest in this figure are the strong differences between the dynamic height current components and between the two sets of ^-hour averages that bracket the time interval when the dynamic height components were computed. 25 The results of these experiments, notably the results of l8 and 20 August, raise strong questions about the extrapolation of current data from any given single current meter station. There are often serious discrepancies between the ^-hour averages in corresponding periods between the two current meter sta- tions that are only one-half mile apart. Further, the drogues move near the current meter stations and show net transports over comparable ^-hour periods that are at variance with the average currents computed at the stations. Finally, the sequence of quarter-period current meter averages shown for the cases of the frontal passages on l8 and 21 August show that the dynamics of the response of currents to winds are not well understood. Moreover, in the period preceding the front of 21 August when the winds were from the southerly direction, the general tendency of the current meter data was to show an overall southward drift. 2k References AYERS, J. C. and ROGER BACHMNN. 1957- Simplified computations for the dynamic height method of current determination in lakes. Limnol. Oceanog. 2:155-157, STERN, MELVIN E. 1965* Interaction of a uniform wind stress with a geostrophic vortex. Deep-Sea Research 12:555-367. 25 III. DETERMINATION OF EDDY VISCOSITY AND EDDY DIFFUSIVITY IN LAKE MICHIGAN Joseph C. Huang INTRODUCTION Theoretically, a flow field of incompressible Newtonian fluid can be solved by the equation of conservation of mass and by the Navier-Stokes equa- tion of motion, satisfying the respective boundary and/ or initial conditions. When the flow is viscous and laminar, a term appears in the equations of motion that is related to shear stress which, according to Stokes' law of viscosity, is du dy where the subscript "l" denotes ''laminar, " y is normal to the direction of the flow velocity, u, and ij, is the molecular viscosity (which is a property of the fluid depending only on temperature under ordinary conditions). In compressible laminar flow, when temperature variations are involved, the heat flux (following Fourier's law of thermal condition) is shown as ^ dT where k is the molecular thermal conductivity, which is also a property of the fluid and usually depends only on temperature. On the other hand, if the flow is turbulent, a kind of motion with random irregular fluctuations is superimposed on the main flow. Since its fluctuations are irregular, it is impossible to describe the motion in all details as a func- tion of time and space. However, one can study its average behavior using the statistical methods pointed out by Hinze (1959). Therefore, in describing a turbulent flow in mathematical terms, it is convenient to separate it into a mean motion and a fluctuation, such as u == u + u T = T + T' 26 where the "bar" indicates mean quantities and the "prime" indicates the fluc- tuation from the mean. The presence of random irregular fluctuations of eddy motion of the flow appear as apparently increased viscosity or diffusivity (Defant I961) . In tur- bulent flow, therefore, additional Reynolds stresses must also be considered in the general governing equations of the motion. Analogous to the coefficients of molecular viscosity, Boussinesq (1877) introduced a mixing coefficient, A^, for the Reynold stress in turbulence ex- pressed in terms of the mean velocity gradients in the flow field, ^ _ du t " " ^ u^v' " t dy where the subscript t denotes turbulent flow. In a similar definition to that for molecular diffusivity, Schlichting (1968) postulated that the turbulent heat flux is t P q qy The exchange coefficients of turbulent momentum and turbulent heat flux, A^ and A^, both have the dimension of a viscosity, \i (g/cm/sec). If the turbulence is isotropic (i.e., the relation between stresses and strains is independent of axis) then A. corresponds to the molecular viscosity \i, A corresponds to the molecular heat diffusivity (k/c ); and e = At/p, and the eddy kinematic viscosity corresponds to molecular kinematic viscosity V = [x/p* Turbulence shear stress is expressed as du T = p € — . t dy The eddy kinematic viscosity, e can be expected to occur as a constant only if the turbulence field is homogenous. Generally speaking, the eddy viscosity is 10^ to 10^ greater in magnitude than molecular viscosity, except in the vis- cous sublayer very near the solid boundary. The latter may be neglected in many cases to a good approximation. As pointed out by Lamb (1932), in describing the wind produced current, \i is replaced by A. . By the techniques mentioned above, it is possible to attack geophysical fluid problems which are usually turbulent in nature. In this way we can probably predict or evaluate the time-mean values of these complicated turbulences whose complete theoretical formulations have so far proved impossible (Lamb I952). 27 Most of the geofluid modelings are treated in this way. For example^ previous work by Sverdrup (19^7)^ Stommel (19^8)^ Munk (1950)^ Munk and Carrier (1950)^ Bryan (I965), and Birchfield (I967) all follow this kind of approach. Nevertheless not all the turbulence is homogeneous and isotropic. The eddy viscosity usually does have directional preference, Eddy viscosity be- comes effective only if there is some flow in the fluids and its value depends on the velocity distributions and on the characteristic length of the flow field. As shown in Boussinesq^s hypothesis and von Karman*s similarity theory t dy 'dy where L is the Prandtl's mixing length. Therefore A^ or e is not a property of the fluid itself^ whereas the molecular viscosity^ \i^ is. It is still extremely useful to introduce the eddy viscosity in relating the Reynolds stress to the mean velocity gradients of the flow. Though A, (or g) is not a physical con- stant^ but varies from one part of the flow to another, estimation of the magni- tude of the average eddy viscosity is still the key factor to the success of geophysical fluid modeling. In Lake Michigan, though the mean velocity is in general less than l/2 m/s (Ayers^ et al. 1958)^, the flows are turbulent in nature. Therefore, it is important to estimate the order of magnitude of the eddy viscosity and the eddy diffusivity in order to facilitate geophysical fluid modeling of Lake Michigan. THE MAGNITUDE OF EDDY VISCOSITY AND EDDY DIFFUSIVITY IN LAKE MICHIGAN In order to evaluate the magnitude of the eddy diffusivity and eddy vis- cosity associated with the thermal current structure of the late spring/early summer circulation in the Lake Michigan model, a series of experiments were carried out to verify order- of -magnitude values for the eddy diffusivity and eddy viscosity to be used with further theoretical investigations. An Estimate of Eddy Viscosity by Reynolds Stress and Vertical Velocity Gradients In August 1967^ a pair of current meter buoy stations were set 8 miles off- shore from Sheboygan, Wis. (see p. 12 of this report for a description of the buoy stations and their locations) » During the period of the current measurement, the research vessel INLAND SEAS was operating around the stations (during the daytime) whenever weather permitted. Wind, as well as other weather data, was automatically recorded on a digital meteorological tape recorder. The anemometer and vane that feed into the recorder are on the mast of the ISLAND SEAS ik m above the water surface. The Captain ^s log and hourly meteorological observations recorded wind from a 28 lower vane and anemometer that are about 11 m above the surface. Using the analyzed recorded mean wind versus height plotted on a log scale with the mean velocity zero at about 1 cm height (Floyd C. Elder, personal communication) the 10-m level wind was interpolated with reference to the Captain *s log. The 10-m daily winds are shown in Table J-l. The vertical velocity gradients in the direction of surface current were deduced from two current meters separated vertically by 10 m. According to Deacon (I962), the Reynold stress is = A ^ t dy • — 3 Then using p^ = 1 g/cm^ and p^^ = I.25 x 10 g/cm^ the eddy viscosity can be calculated. The data^ wind, wave conditions, absolute vertical velocity gra- dient, shear stress, and calculated eddy viscosity are all shown in Table 5-I. The values of Reynold stresses obtained here using Deacon's formula are in good agreement with the mean value of stress by Elder and Soo (1967) at The Univer- sity of Michigan Research Tower on the east side of the lake. An Estimate of the Upper Bound Value of Eddy Viscosity by the Result from Statistic Theory It is quite well known (e.g., Ayers, et al. I958) that a large clockwise eddy exists in southern Lake Michigan between Grand Haven and Michigan City. This large eddy has the length scale of a little less than half the dimension of the width of the lake. Bowden (196^1) has shown the effective eddy viscosity is proportional to the V5 power of the length scale, L, of the particular eddy under consideration. In Lake Michigan, let this largest eddy have average length scale of L = ^0 km or L = i+ x 10^ cm. Then, according to Richardson's formula, the eddy viscosity will be of the order lo'^ cm^/sec which is the upper bound of the numerical values of the eddy viscosity in Lake Michigan. An Estimation of Vertical Eddy Viscosity from Current Measurement Near the Bottom From 11 May to 11 October ±9^7 y a tripod- supported pendulum current meter was installed about 25O m offshore, T-l/2 miles south of Benton Harbor in a depth of 4.5 m. The pendulum was suspended only 20.3 cm above the sandy bottom. The monthly mean current is shown in Table 5-2. Using Lesser 's (195I) roughness length for sand bottom, z^ = .I5 cm, the friction velocity near the bottom can be determined through z + z \ = V[in(-^)]"' o 29 (D H •H B 00 -P q5 VD ON ■P S, (M I ^ O •H ■P cd ■P w o -P I +:> (L> O o o ■P •H Ul O O m > cd t>. o rQ CD CO >? Id Ch H O H I pq >. +:> •H w ^-> o o O 0) w w o •H tsl \>. (M > <; CM^ VO g >= o T^ v_. ^ w H Ul * O W CD W CD -P Ph a:! CO >^ 0. O •H -P O (D CD ■P cd cvi I VD CJ\ CM \\ c^j H i bO O H CM CO CM CM H 00 H CM CM CM I O H VD H CO CM 's 1 O H X XI 00 00 H LTN CM CM H CQ o CM t-- H CM r^ CO CM CO I O H H CM CO m Ul CQ »- CO (D ir\ CD hTN CD g LTN CM CD > 1 !> 1 t> H 1 1> 1 fe cd - cd - Cd cd - Cd o.- Cd > N^ ^ H ^ o K^ > H ^ o\ CO bO CM X (D • a H o cd > "So CD CO bD cd b ^ H CD 1> X! cd K^ CM 30 where u is the mean velocity, z is the height from the bottom during measure- ment, and Kq = O.Ul, The monthly mean friction velocity for z = 20.0 cm is shown in Table 5-2. The neutral stability condition generally persists near the bottom. Then the vertical eddy viscosity can be estimated as shown in Table 5-2 (Bowden I96J4) . TABLE 5-2. Estimated vertical eddy viscosity by current measure- ment near bottom. Vertical eddy Mean current Friction velocity viscosity (cm/sec) (cm/sec) (Az = k u z cm^/sec) Month o^ May- 59.5 June 51.7 July 1+7.2 Aug. 58.5 Sept. 56.1 Oct. 69.8 k.6k 2.krj 5.68 I+.56 U.58 5.UU UO 20.7 50.9 38.5 56.0 1+5.7 Average value 33 cnf /sec An Approximation of the Vertical Eddy Viscosity by Surface Wind Wind data from both the two- current-meter-buoy operation off Sheboygan and the wind and current station outside of the Benton Harbor beach are used for approximating the vertical eddy viscosity in Lake Michigan. In the two- current-meter-buoy operation off Sheboygan, wind data were measured aboard the R/V inland seas. At the current meter station outside of the Benton Harbor beach, the wind was measured Just above the water. From Sverdrup, et al. (I96I+), depending on whether the magnitude of wind was greater or less than 6 m/s, the vertical eddy viscosity was approximated as shown in Table 3-3. Notice that the vertical eddy viscosity is smaller when the wind is averaged for a longer period of time . 31 TABLE 3-3. Estimated vertical viscosity by wind speed. Vertical eddy Time Wind (m/s) viscosity, A^ cm^/sec Remarks IT Aug. 1967 9.5 18 Aug. 1967 9.8 20 Aug. 1967 5.7 21 Aug. 1967 5.^ 25 Aug. 1967 6.2 2k Aug. 1967 h.6 581 165 100 Wind averaged over 6 hours Average value 206 cm^/sec May 1967 3.5 June 1967 2.5 July 1967 5.1 Aug. 1967 5.2 Sept. 1967 2.8 Oct. 1967 k.l Nov. 1967 k.k i;^ Monthly mean wind 16 50 55 22 TO 87 Average value ^3 cm^/sec A Measurement of the Eddy Diffusivity by the Dispersion of Marked Particles In a series of experiments conducted in Lake Michigan, the fluorescent dye rhodamine B was used as a tracer for diffusivity measurements. Rhodamine B has peak ultraviolet fluorescence wavelength of 58O millimicrons. In May I968, a fluorometer with only one sample cell was used. Because in this season the phytoplankton was so abundant, and since they have an ultraviolet fluorescence frequency Just a little lower than that of rhodamine B, the fluorescence of the peak concentration of rhodamine B dye was not distinguishable from the background after an hour or so. In all subsequent experiments, a comparative measuring fluorometer with two cells was used. The rhodamine B mixed with natural lake 32 water flows through the dye sample cell, while the standard cell contains natural lake water. The relative difference in fluorescence between the two cells is measured in order to obtain reliable data in the possible presence of any inten- sity variation of ultraviolet light lamp,^ systematic changes in the characteris- tics of photocell, and differences in the amount of naturally occurring fluorescent substances in the water (Noble and Ayers I961) • Both instantaneous source and continuous source measurements have been made in I968. Instantaneous Sources . Three runs of instantaneous sources of marked particles for diffusivity studies were conducted. These were 1 and 8 miles off Grand Haven in May, and 5 miles off Milwaukee in June. Rhodamine B dye was diluted with methanol to have a specific gravity very close to the lake water. Generally, one part of stock solution of rhodamine B is Uo^ acetic acid solu- tion was mixed with 5-1/^ parts of methanol. The dye mixture was gently poured out on the surface lake water at the upwind side of the boat in order to allow the boat to drift away with least interference to the dye patch. The first and second runs were investigated aboard the r/v MYSIS and the third run aboard the R/V INIAKD seas. During the first 20 minutes after the source had been instan- taneously releaseid, the diameters of the expanding dye patch were estimated by eye with the ship's length as reference. At 20-minute intervals, the boat was used to coast through the dye patch longitudinally and laterally at the lowest possible speed. For each transect, fixes were obtained from shore objects when entering and leaving the patch, along with the time required to traverse through the patch. Continuously running back and forth across the dye patch, the con- centration of the dye patch was continuously recorded from the fluorometer on an AZAR recorder. The concentration distribution of a typical pass through the patch is shown in Figure J"!- Most generally, the dye concentration distribution showed a form of Gaussian distribution. The diameter of the patch, longitudinally or laterally, is properly approximated as k times the standard deviation of the patch from its mean position. Batchelor (1955) shows that (assuming the fluid is in a uniform, homogeneous, and steady, mean velocity field) at the initial stage of diffusion process, the effective diffusivity (which is proportional to the rate of change of the square of standard deviation with time) is proportional to the first power of the elapsed time; in the intermediate phase, it grows as a V5 power of the patch size; and at the final stage, as the standard deviation grows as l/2 power of the elapsed time, the diffusivity should asymptotically approach a constant. The general information and the value of effective diffu- sivity are shown in Table J-U. The growth of dye patch versus time is shown in Figures 5-2, and 5-5. Continuous Source . Four runs of continuous (in time) point sources of marked particles were conducted near the meteorological research tower 1 mile off the east shore in Lake Michigan. Ten gallons of rhodamine B and methanol mixture were contained in two 5-gallon carboys connected together by plastic tubings as shown in Figure 5-^. The continuous point source was released from the apparatus through a rubber tubing which had its outlet about 20 cm above the water surface. The rate of discharge of the dye was controlled by a screw clamp on the tubing,, and was usually allowed to flow at a rate of 5 gallons/hour, or approximately h g/sec. In a homogeneous fluid with a steady uniform current, the continuous dye source was observed as a slender plume extending downstream. 55 r «' 4' Distance (in unit scale) 6^UT Figure 5-1- A typical distribution of dye concentration across a dye patch. TABLE 5-i^■. General information and effective eddy diffusivity ©f dye patch studies. Run Run 1 Run 2 Run 5 Date 1 May 1968 2 May I968 27 June 1968 Location 1 mile offshore 8 miles offshore k miles offshore near Grand Haven near Grand Haven near Milwaukee Depth 20 ft ko ft 65 ft Duration of run 1150 - 1?10 EST 1500 - IU50 EST 0910 - 12if0 EST Current velocity cm/sec 11.6 5.5 7.6 Wind Direction m SW SSW Speed (kts) 15 15-20 5 Lake condition 2-5 ft waves 1 ft wave 2 ft waves Effective eddy dif- fusivity (cm^/sec) Lateral Not available lli^O 7^5 Longitudinal Kot available ^131 ^-019 Remark Very poor visi- bility 34 l/>00< 600 300 too 60 30 ,T^ .rife f / y / / / / / (a) T6 '1.2 • ■3.0' ' ■ ' "6 Time (minutes) 1,000- 600 30 0-J I o too 60- / y/z / / / / / X^ / / X/ (b) n:! ^ ' *3^' ' ' ' '6 Hi ' ' '30 60 ^20"^' '- Time (minutes) Figure 3-2. Lateral growth of the dye patch with time, (a) 2 May I968. (b) 27 June 1968 55 I,000- /* 600 / 300 / / / 4J 0) / w / ^ X 100 / / / SO- SO 10- (a) 600- 100- Time (minutci} / X / / X / / / / / / / / x/ / / (b) .- Ti—" — ■ "s-d ■ ■ ' V. ■ "« ■ 'so' ' ■ ' "60 "So"" — ' Time (minutes) Figure 5-5, Longitudinal growth of the dye patch with time, (a) 2 May 1968. (b) 2J June I968. 56 Air Intake- Screw Clapper* Figiire 5-^. Continuous source apparatus. The direction of the dye plume was reported by the observer using a theodo- lite mounted on the meteorological research tower 57 ft above the water surface. One or two buoys were dropped as markers on the edges of the dye plume at right angles to the direction of flow. The appropriate distance of the marker from the source depended on the mean velocity of the current which was estimated from the drogue study set earlier. The r/v MYSIS crossed the dye plume back and forth at 10-minute intervals at the lowest speed she could make (about k.^ mile/hour). As soon as the dye plume passed through the markers, samples of water were taken through a constant- flow fluorometer that fed into an AZAE recorder. Communication between the tower and the boat was maintained through a VHF transceiver. Whenever the boat entered and left the dye plume a "MARK" signal was transmitted to the theodolite observer. The observer obtained the azimuths of both edges of the plume and the elevation from the tower to the boat. The horizontal distance between the boat and the tower was obtained from the elevation angle as well as by fixes from land objects. The angle subtended by the edges of the plume is the difference of the two azimuths. The arc length of the plume at this distance from the tower was arbitrarily taken as four times the lateral standard deviation unit. At each return crossing, the vertical diffusion was measured by lowering the fluorometer intake at half -meter inter- vals until there was no indication of the dye. Half of the maximum vertical displacement was taken as the vertical standard deviation unit. 57 At the last few minutes before the exhaust of the point source^ continuous crossings were conducted at various distances along the plume from the tower to approximately I-I/2 miles downstream. These crossings measured the asympto- tically steady diffusion of the plume. These continuous point source diffusion studies were carried out on 50 July and 1-2 August I967. On 50 July, sampling was discontinued about 2 hours after the experiment was initiated, because of high seas. The dye plume was contin- uously sampled for 5 hours on 1 August, and during the sampling period the cur- rent shifted its direction of flow about I5 degrees. Two sampling runs were conducted on 2 August under a very steady current condition. The general infor- mation for each run and the value of diffusivity are shown in Table 5-5. The lateral growth of dye plume versus time at a constant distance from source is shown in Figure 5-5. The lateral growth versus distance under asymptotically steady condition is shown in Figure 5-6. A Measurement of the Horizontal Diffusion by Drogue Studies A drogue is made of k pieces of light sheet metal (e.g., aluminum alloy) approximately 1 x 2 m in dimension, linked together in the shape of a cross and suspended by chains in the water by a float which supports a radar antenna above the surface (similar to Farlow I965). It is frequently used for measuring the mean current of the flowing layer at the drogue depth. The equivalent radius of the drogue is about 1 m and, based on its geometry and the corresponding Reynold numbers (10^ to 10^ in the Great Lakes), its time constant is estimated to be 100 sec.(Okubo and Farlow I967) . Though the small-scale turbulence may be averaged out, the drogue is adequate for measurement of large-scale turbulent diffusion for eddies larger than 10 m in scale which represent the energy- containing eddy scales in the Great Lakes (Csanady I965; Okubo and Farlow 1967) . On 27 June I968 a hexagonal pattern of surface drogues, with the center one very close to a fixed buoy, was set about 5 miles offshore near Milwaukee in Lake Michigan. The original distance between each pair of drogues was l/2 mile. After setting the pattern, the r/v INLAND SEAS continuously fixed each drogue for the following 6-I/2 hours. Similar to the result mentioned by Okubo and Farlow (1967), the growth of the pattern size with time was not clearly noticed. The diffusion constant was deduced from the second moments about the mean of the haxagonal pattern. At any particular time, the instan- taneous positions of all drogues were interpolated from their respective pre- vious fixes. The general information on the experimental run and the value of eddy diffusivity are shown in Table 5-6. COMPARISON OF THE VALUES OF EDDY VISCOSITY AND EDDY DIFFUSIVITY WITH OTHER STUDIES In Lake Michigan, there is little data about eddy viscosity or eddy dif- fusivity. Some figures from the ocean or other lakes are available as reference « 58 CD O U g w § •r-l -P o o O >:. -P •H t> •H W ;:^ Ch Ch •H <D CD > -P O CD Ch Ch CU cd s^ O o • Ch >5 -P m cd 0) es ft pq CD •s EH iH CD CO ^ Sh H CX) ^ •H i) MD Ch •H -d- ON Ch cd •n ^ LTN H § § •H H • CD ^ bO S 1 -d- tiD H -P •H ?H r^ CD Ch j;^ < S C^J W M H J2J CD 5 cd CM « CM H CD CD S h^ \s\ H EH M fH pq XI c6 a Cd CX) w Hi •H W LTN VD Ch •H ^ C^ Ch cd •\ ^ On H g g » CD X bD 1^ 1 hr\ tiD H CD -P ^ 'd U M CD Ch 5;;^ < n cd OD ^ jj CD 5 Cd CX) p:^ CVJ H CD a? S H^ Lr\ EH CO X S W ^ cd a cd 00 Ul h-:) •H bD 1^ VD Ch •H hf^ o^ Ch cd •\ ^ hCN H Cd H • S bD :^ 1 CM bO H CD -P jii •H ^ r^ <D Ch \r\ a «^ S cd w ^ H j:^ CD 5 •^ pc; H H CD s h^ ir\ EH CO >:! a H 00 ^ cd a cd VD w i-q •H bO ON Ch •H CM H Ch cd •\ ^ CM 0. a H !>5 cd •H H CD S bO S 1 H H CD P> h) ^ M CD Ch f^ ^ cd CO cM a jj CD jj Cd r— P^ rr\ H 0. I^ i-:| Lf^ r— a •H -P CD Cd P> Cd P i-q ■p Pi CD Ch o o -p cd Pi 00 CM CO CM 00 CM l:^ 00 ^ 00 -p o ^ O TJ 1:1 CD •H P CD CD ft CO CD > Cd ■p Ch CM I cd -p Ch I CM CO CD !> cd 15 -p Ch \r\ w 1 CO H hr\ o •r-l -P •H t;} a o o CD cd P CM CO (D !> cd ^ -p Ch H CM I>- 1 hCN rTN CO Lr\ X^ CD bD Q) g ;:s H •H rQ -P 0. cd CD H bD •H 0. bD Cd •H a % CO H -P cd -P CD ^ X s::! CD P> CD -H ^--N > > O •H 'H p> CO CO O p( \. CD Ch c\j_, Ch Ch g Ch -H o H Id ^^ CD p:; 39 1.000 900 800 700 600. 500- u (U 4-1 400- >> 300- 200 100 / / / / 60 120 (a) '240 ' ^363 ' ■ '600 400 300- 200- 100- 90 80-1 70 60^ 50 18 / ^ .X' '30 ■ '42' • '60 Time (minutes) (b) '120 ^I80" Figure 5-5. Lateral growth of dye plume at a constant distance from the source, (a) 1 August I968. (b) 2 August I968. ko \ \ O v CD O 1 I r (53aj) aouBiSTQ pasnjjTQ T"35B1 \, \. \. is \. \. \ \ \ 1 i i § i i" « rt N M - = (53aj) 33UB5STa pasnjJTQ TW35B1 o o o cd -P m -P •H ^ • B 00 3 vo H ON ft H *i o <ci ^ CM ■P o ^ H ^^ cd 0) CO +3 vo cd ON h^ H -P • CQ I <t) ^ H a §,. •H (d s s § s ♦ •" lO «• (3333) souBqsja P3snjjTa TBa35Bi 41 TABLE 5-6. General information and eddy diffusivity of drogues study in Lake Michigan. Date: Location: Depth of water: Depth of drogues: Duration of run: Current velocity: Wind: Direction Speed Lake condition: Mean initial separation: Final deviation from mean: Effective eddy viscosity: 27 June 1968 5 miles off Milwaukee 65 ft 5 ft I5OT - 1850 EST 7.6 cm/sec N 6 kts 1/2- 1 ft waves 2578 ft 3UOO ft 5A X : 10* cm^/sec In modeling the Gulf Stream, Munk (1950) estimated the lateral eddy viscosity of 5 X 10'^ cm^/sec. Stommel (1955) estimated the eddy viscosity of Gulf Stream near the Florida Straits to be 2 x 10^ cm^/sec. Defant (I96I) determined the exchange coefficient of lateral large-scale turbulence in ocean phenomena as 5 x 10'^ cm^/sec. In general, the vertical turbulence scale is much smaller in dimensions than that of the horizontal. Therefore the magnitude of vertical eddy viscosity is much smaller. According to Defant (1961), Sverdrup obtained the value of the vertical eddy viscosity of 75-260 cm^/sec in the North Siberian Shelf -re- gion; Fjeldstad estimated the vertical eddy viscosity in the Caspian Sea of 0-22if cm^/sec, and Suda measured the value of 680-7500 cm^/sec for the Kuroshio and of 150-1^^60 cm^/sec in the Japan Sea. For vertical eddy diffusivity, some typical values have been summarized by Defant (I961) as: Mediterranean, k2 cm^/sec; California current, 50-^0 cm^/sec; Caspian Sea, 2-l6 cm^/sec; Equatorial Atlantic Ocean, 520 cm^/sec. Csanady (1965,196^,1966) has measured the eddy diffusivity in Lake Huron and Lake Erie, and found horizontal diffusivity of the order of 5OO-25OO cm^/sec in Lake Huron k2 and, 1000-5000 cm^/sec in Lake Erie. These values are much greater than verti- cal eddy diffusivity, which is only of the order 0.1-10 cm^/sec. Okubo and Farlow (1967) measured the horizontal diffusivity for large eddies in Lake Michigan and Lake Erie using drogues and found the eddy diffusivity of 5 x 10 to 6 X lO''^ cm^/sec. According to our estimation, the vertical eddy viscosity is on the order of 55 to 2 X 10^ cm^/sec in Lake Michigan. If the vertical eddy viscosity is usually two orders smaller than the horizontal value (Defant 196l)jj then the horizontal eddy viscosity should be in the range of 3.5 x 10^ to 2 x 10^ cm^/sec, which is within the upper limit of lO"^ cm^/sec as we have estimated. The lateral eddy diffusivity we found using marked particles is 750-1200 cm^/sec for small eddies and ^.k x 10^ cm^/sec for larger eddies. The effective longitudinal eddy diffusivity for small eddies is in the order of ^ x 10^ cm^/sec. All the data we used appear quite consistent and in good agreement with that of others. CONCLUSION In modelling a geofluid problem^ the eddy viscosity and the eddy diffu- sivity are of critical importance in order to predict the natural current or wave phenomena with a similarity to laminar flow. In Lake Michigan, though the mean current velocity is small in general, the flow field is turbulent in nature. In solving or explaining the flow pattern of the mean lake current, it is possible to use the governing equations of laminar flow with eddy viscosity and eddy diffusivity in place of the molecular viscosity and the molecular dif- fusivity* A series of experiments have been conducted in order to estimate these two values. The data presented are quite consistent and in good agree- ment with the data reported by other investigators. The vertical eddy viscosity in Lake Michigan is in the range of 1 to 10^ with a mean 10 cm^/sec. The horizontal viscosity is in the range of 10^ to 10"^ with a mean value of 10^ cm^/sec. The eddy diffusivity may reach the same magnitude as the viscosity but it is in general smaller. A typical mean value for vertical eddy diffu- sivity is 5 cm^/sec^ and for the horizontal eddy diffusivity 10^ cm^/sec. Though few experiments are far from enough to state the characteristic values of turbulent phenomena, the data presented together with values reported by others can be used as good reference values for Lake Michigan. h3 References AYERS, J. C, D. C. CHAKDLER, G. H. LAUFF, and C. F. POWERS. I958. Currents and water masses of Lake Michigan. University of Michigan^ Great Lakes Res, Div. Pub. No. 5, 169 p. BATCHELOR^. G. K. 1955- The theory of homogeneous turbulence. Cambridge Uni- versity Press, Cambridge, Mass. 197 P* BIRCHFIELD, G. E. 1967. Wind-driven currents in a long rotating channel. Tellus 19(2):2i|-5-2i|-9. BOUSSINESQ, J. 1877. Theorie de l^ecoulement tourbillant. Mem. Prec. Acad. Sci., XXII, I4-6, Paris. BOWDEN, K. F. 1964. Turbulence. P. II-50 in Oceanogr. Mar. Biol. Ann. Rev., Vol. 2. H. Barnes, ed. George Allen and Unwin., Ltd., London. BYRAN, K. 1965. A numerical investigation of a nonlinear model of wind-driven ocean. J. Atmos. Sci., 20:59^-606. CSANADY, G. T. I963. Turbulent diffusion on Lake Huron. Jour. Fluid Mech., 17:560-585. . 1966, Accelerated diffusion in the skewed shear flow of lake currents. J. Geophys. Res., 7I: 4ll-J^-20. DEACON, E. L. 1962. Aerodynamic roughness of the sea. Jour, of Geographic Res., 67:5157-5172. DEFANT, A. I961. Physical oceanography. The MacMillan Co., New York, Vol. 1, 729 p. ELDER, F. and H. K. SOO. I967. An investigation of atmospheric turbulent transfer processes over water. University of Michigan, Great Lakes Res. Div. Spec. Rep. No. 29, k6 p. FARLOW, J. S. 1965. A field technique used for horizontal diffusion studies in Lake Michigan and Lake Erie. University of Michigan, Great Lakes Res* Div. Pub. No. 15, pp. 299-505. HINZE, J. L. 1959. Turbulence. McGraw-Hill Book Co., New York, 586 p. LAMB, H. 1952. Hydrodynamics. 6th ed., Dover Pubs., New York, 758 p. LESSER, R. M. 1951* Some observations of the velocity profile near the sea floor. Trans. Am. Geographic Un., Vol. 52, pp. 207-211. kk MIMK, W, H. 1950. On the wind-driven ocean circulation, J. Meteor,^ 7(2): 79-95- and G. F. CARRIER. 1950. The wind- driven circulation in ocean basins of various shapes. Tellus 2(5) : I58-I67. NOBLE, V. E. and J. C. AYERS. I96I. A portable photocell fluorometer for dilu- tion measurements in natural waters. Limnol. Oceanog. 6:U57-^6l* OKUBO, A. and J. S. FARLOW. 19^7. Analysis of some Great Lakes drogue studies. Proc. lOth Conf. on Great Lakes Res., Internat. Assoc, for Great Lakes Res., pp. 299-508. SCHLICHTING, H. I968. Boundary- layer theory. 6th ed. Translated by J. Kestin. McGraw-Hill Book Co., New York, jkj p. STOMMEL, H. 19^8. The westward intensification of wind-driven ocean currents. Trans. Am. Geophys. Union 29(2) : 202-206, SVERDRUP, H. U. 19^7. Wind-driven currents in a baroclinic ocean; with appli- cation to the equatorial currents of the eastern Pacific. Proc. Nat. Acad. Sci. 55(11): 518-526. ^5 IV. MJMERICAL EVALUATION OF STERN ^S CIRCULATION MODEL James H. Saylor DEVELOPMENT OF STERN^S MODEL In 1965^ Stern presented a model for the interaction of a uniform wind stress with a geostrophic vortex. The results of these calculations showed that a geostrophic vortex subjected to a uniform wind stress would move, as an entity, with an Ekman drift, provided that the thermocline was "soft" and that internal waves could be generated at the center of the vortex. Stern* s results have been evaluated in terms of wind stresses and vortex dimensions that might be expected to occur in Lake Michigan (Noble 196?). The results predict a very broad range of internal wave periods that are within the limits of values measured for Lake Michigan (Mortimer 1965). The starting point for the theoretical model is the assumption of a hydro- static, two-layer model; upper layer of density p and lower layer of density p^; no friction; and the fluid incompressible. A A Dv J ^ VP A J1 — - + fkxv+ — = • -B'^*^ Dt p (1) where H = depth of upper layer. and p - p o D d ^ ^ A A A A A U = V + Wk = Ui + vj + Wk A V-U = Let f/h = (S/Sz) G where 9 is permitted to be a lateral wind stress only. ^In par- ^ / \ A ^ ticular, take (x,y,z = o,t) equal to the surface wind stress, t. Let U = Uq + U-j^, where 11^ is taken as the geostrophic component of (l). k6 A ^ h U = V + KW o o o A / A u = V + m^ b b b Thus A D V o o Dt 'op A /I VP / + fk X V + — = - g'k (2) o where A V-U =0 , o A A and whence^ U-j^ is readily defined as the component of U due to the surface wind stress, T. Call this part the "Ekman drift." Boundary conditions are established by requiring no displacement of the free surface^ i,e,^ ^'U^ = - k^^U^ = W^(x^y^z = O)^ and by requiring at the fluid interface that the "Ekman drift" is limited to the upper layers^ and that the wind current joins the geostrophic component smoothly Just above the interface, ^A U^(x,y,z = -H,t) =-^(x,y,z = -H,t) = , which implies 5. ..^ ^ ^^ n 9(x,y,z = -H,t) = ^ = ^ = A/ s A 9(x,y,z = 0,t) = T Subtracting the geostrophic component, equation (2), from the momentum equation, equation (l), gives: A which can be written as: b^AAAAA Aa^G ,, -T— + V • Vv + V -Vv + V -Vv + f k x V = -:r- . ( 5 M at boobbb bdz ^^ ^ hi Taking the vertical component of the curl of equation (5')^ denoting: if s k.V X U = k-V X V ^o o o A 4 4 A Stern obtains, where: ^ A /I /I A 2^ = k-V X U = k*V X V ^b b b (f + <f )-^— + -^(k-V X 0) = v^-ViT + Z T. (4) o' az 3z' b ^o i ^ ^ 1=1 dz dx o b dz dry' o b o D o b Ts Sz 5x dz ^ ' a <^ ^ .. ^f Ts = (W + W^)^— + W^ . o b dz b dz Integrating equation (h) between z = -H and the free surface: /^ (f + ^ )^ dz + /^ |-(k-V X 0)dz = /^ v^.V^'dz + /^ Z T.dz (5) -H ^ ^o^Sz -H Sz^ ^ -H b o -H . , i ^^^ 1=1 which becomes: •(f + W )W^ - k-V X T + f V -V? dz = - / Z T.dz . (6) o -H b o -H . , 1 1=1 i+8 stern compares the magnitude of the terms in equation (6) and shows that if the Rossby numbers based on both the geostrophic and "Ekman" currents are small, i.e. y D ^ O _ the / T. terms are of negligible magnitude. Defining the wind-driven volume transport as; /^ v^dz = M = - - k X T (which is the Ekman transport in the conventional sense). (f + ^ )W^ - k-V X T + M-V^ = (6') where A a h ^ ^^ k X T-V6 To include the effect of f varying with latitude, use the beta-plane approxi- mation, which simply replaces the V<^ term with V(f + ^ ). So, the general so- lution is written: which may be expressed as: «* = -'-|ff|^ (T") where For the vorticity of the geostrophic component in the upper layer we have the relation: ^9 Integrating this expression between -H and the surface: ,0d,. .,, ,0 ^ Ih H(' ^ i'- = /^ S^ - ' A A dt which may be written as: d(f + ^ ) dt H , f rr ^ X T f dH f + £' ~ H dt (9) Equation (9) is the generalized equation of motion for the geostrophic component in the mixed "Ekman layer." For a constant, uniform wind stress (v x r = 0) and for a limited area in which f ^ constant, equation (?') reduces to: (k X T)*Vf ^^=-77— F^ (10) ^^-o o and the kinematic equation, equation (9)i becomes: H dt " H ( f + f' )^ dt . (11) Stern's analysis with the necessary conditions and restrictions, describes the interaction of a localized, geostrophically-balanced flow (i.e., within the confines of the simple, two-layer model envisaged) with a wind-driven "Ekman" current. The vertical velocity (W^) generated by this interaction is propor- tional to the gradient of the geostrophic vorticity (V^), a quantity which is conserved in the upper layer, and to the wind stress, t. In this manner, it is seen how the vertical velocity transfers momentum across the interface to the lower layer, and how there is a tendency for the vorticity of the water column to increase due to the advection of neighboring vorticity by the mean Ekman drift. 50 APPLICATION OF STERN'S MODEL TO LAKE MICHIGAN If the wind stress (t) and the coriolis parameter (f) are independent of position, the last term of equation (ll) reduces to: ^(kXT)-V^^ ^k.(V^^XT) k.Vx(^^fT) a k.v X (— ) which is equivalent to a force per unit mass of - ^^r/Ef acting on the upper layer. Consider a two-layer, model (Fig. k-l) with Ekman drift at right angles to T and directed along the x - axis, i.e.. T = positive constant y T = 0. X ^ ^_^ ^ ^ ^ ^ -^ ^ ^ ^ - • ^ ^ X^ Vorticity <, I E^ /> .y -\ + \ /// ^ /////// //> ^ ' Vorticity li a. ^ H^ '• ^ ' 5- -H, - H, Figure i|.-l. Schematic two- layer model. Let the positive displacement of the interface at z = -Hi be small and given "by Ti(x,y, t). The vorticity equation for the lower layer is: d4 _ ^ ^ dt az • Integrating between z = -H^ - Hi and z = -Hi + r| ,-Hi+Ti d^ _ pHi+n ^ ^ ^^ iHi-H2 dt ^^ iHi-He ^ hz ^^ ' dt Ha' 'z = -Hi+Ti , 51 ^ ^z = -Hi + n dt Integration gives: ^*#- = ^*it" ' t=0 t=0 ^^(t) = Ul^llll . |^^(^^y^t=0) . (12) The vorticity equation for the upper layer is: 'i^^i'^^^^t'" where U = r/fHi = mean "Ekman" drift. (15) The interaction of an Ekman drift current with a geostrophic vortex can be described in terms of a vertical suction velocity W^, which is proportional to the gradient of the geostrophic vorticity and to the wind stress. This interac- tion leads to an increase in the total vorticity of the water column with time because of the transfer of momentum across the interface. It was shown that if T and f were constant, the result of the interaction can be described as advec- tion of vorticity at the rate of mean Ekman drift. Equation (I5) expresses this fact. If the geostrophic vortex is centered at x = y = at time t = 0, at which time an Ekman current is established, equation (I5) states that part of the vorticity will stay at x = y = (represented by the interface disturbance), while part is advected with the drift current. An expression for the vorticity difference between the upper and lower layers comes from Mar gules * formula *ri =<& - f^vit) (Ik) where r^ -^ ^ ^ -^""W From (12), (15) and (lU): /l. + ui-)(^ .^Ivin) + — -^ = U — ^ ^ ' "^ hi'^-YLz f '^"'^ ' Hi ^ "" H2 ax ' ^"""^^ 52 Separating equation (l4) into steady-state and time- dependent components by- setting Tl(x,y,t) = Ti^(x^y) + 7i^(x,y,t) , for the steady- state component: for the time- dependent component: Assuming a solution of the form: ik(x-Ct)+i% ''t " ° (C = phase speed) then \ = ^ vin^ = -(k- + i^)!^^ equation (I7) becomes: ^^^Hi_^^ + g' p(k^ + -«^)] = U[l H- g' ^(k^ + i^)] , (17.) or, §=[l+|(l-k--^0"" (18) where K = (^)'/^ k ^, = (^)V^ i 55 -4 -1 Take f = 10 sec (an approximate value for Lake Michigan). If k^ » Hg/Hi, C 21 U and short waves are advected at the speed of the mean Ekman drift. If k^ = i^ is small- the phase speed is: (1 + k^ + 4)Hi and a surface vortex of radius ^(g* -3)^/^ will leave the origin in time: ^"2^ f^ ^ ' LuHi(l + 4 + 4) J • To estimate the eddy diameters, ratio of propagation, and internal wave periods involved, use the following parameters as indicated in Figure i|— 2. I - ° a j^ = 20 m ^ = 0.9982 (20° C) H2 - 100 m /tf© s^+A^= 0.9998 (8" C) Figure k-2. Specific parameters of the two- layer model. Let: ^:^ = ^:^ ^ f = 10" sec" :' = (— )g « 1.57 X 10"^ m/sec^ T = 2 X 10' (m/sec)^ Then, *T* U = — - = 10 cm/sec rHi % ^ ~ ^[{1 H- k^ + 4)Hi + hJ eddy diameter = — ( ^g^ )"'"/^ TABLE i+-l. Calculated:eddy diameters, phase velocities, and interval wave periods of upper layer in Lake Michigan. K = h Eddy diameter, km cm/sec Phase speed c, cm/ sec Wave period, days 1 59 10 5.75 12.1 2 1.95 10 5.0 ^.57 5 15 10 6.67 2.25 it- 7.8 10 9.1 0.99 5 ^.8t 10 9.55 0.59 Effect of Density Difference, Ap Keeping other parameters fixed, let Ap = l/2 of former value, i.e., take g' = 0.8 X 10 (cm/sec^). Considering the density change, the result is cal- culated as shown in Table 4-2. TABLE 4-2. Calculated eddy diameters, phase velo- cities, and internal wave periods for upper layer in Lake Michigan. Eddy k = i diameter, U, km cm/ sec Phase speed cm/sec Wave period, days 1 28.1 10 3.75 8.7 5 9A 10 6.67 1.63 8 5.5 10 9.55 O.J+5 55 Effect of T The wind stress determines the speed of the Ekman drift. Within the limits of (v-[^|fL « 1 then^ the only effect of t is to determine the rate of advection* If the parameters of the basic model remain fixed, and t is reduced to l/2 of its former value (1 x 10 (m/sec)^), the solutions are modified as shown in Table k-3. TABLE ^4-5. Calculated eddy diameters, phase velo- cities, and internal wave periods for upper layer with wind stress changed. Eddy ^^ ~ ^^ diameter, km Phase speed Wave U, C, period, cm/sec cm/sec days 1 59 5 1.87 24.2 5 15 5 5.55 i^.5 8 k.Qj 5 U.76 1.18 Effect of Depth of Thermocline Assume the same parameters and overall depth as basic model, but let Hi 10 m, H2 = 110 m, then the results will be as shown in Table h-k. TABLE k-k. Calculated eddy diameters, phase velo- cities, and internal wave periods in Lake Michigan with the depth of thermocline changed. K = h Eddy diameter, km cm/sec Phase speed c, cm/sec Wave period, days 1 5 i^5.4 14. U 20 20 1 4.28 12.6 11. T 1.52 For Hi + H2 = constant, the effect of decreasing Hi is to increase the mean Ekman transport velocity (u), to increase the wavelength associated with 56 each wave number^ and to decrease the wave period (because of increased advec- tion rate) . The range of values for internal wave periods calculated from this model are within those reported by Mortimer (I965) and also within the range of perio- dicities of current speed fluctuations calculated by Verber (1965). However, because of the limitations imposed by the assumptions of a simplified model, these values cannot be considered to be definitive. Stern* s model does provide mechanisms for the generation of phenomena observed in field measurements, such as rotary currents from buoy stations, rotary drogue motions, internal waves, and a transfer of wind energy into the water circulation in a way that might explain the buoy station measurements re- ported in the previous section of this report. It will be necessary to design future field experiments relating the three- dimensional temperature field to the three-dimensional current field as a func- tion of time. These experiments will be time-constiming and expensive, but the numerical results presented in this paper indicate that there is a strong reason to make definitive tests of the Stern model. 57 References MORTIMER^ C. H. I965. Spectra of long surface waves and tides in Lake Michigan and at Green Bay^ Wisconsin. ProCo 8th Confo on Great Lakes ReSo ; University of Michigan^ Great Lakes Res. DiVo Pub, No, 15:50^1-525. NOBLE^ V. E. 19670 Evidences of geostrophically defined circulation in Lake Michigan. Proc. 10th Conf. Great Lakes Res«^ Internat. Assoc, for Great Lakes Res., p. 289-298. STERN, MELVIN E. 19650 Interaction of a uniform wind stress with a geo« stropic vortex. Deep-Sea Research, 12:555""67o VERBER, JAMES L. 1965» Current profiles to depth in Lake Michigan. Proc, 8th Conf. on Great Lakes Res. ; University of Michigan, Great Lakes Res. Div. Pub. Noo 15:56^-571. 58 V- THERMAL CURRENT STRUCTURE IN LAKE MICHIGAIT, ASSOCIATED WITH THE SPRING WARM-UP SEASON— A THEORETICAL STUDY Joseph C. Huang INTRODUCTION The early spring condition of the lake lends itself most readily to mathe- matical analysis. During the winter^ the water mass of the lake becomes nearly isothermal^ and has a three-dimensional circulation pattern reaching to the bot- tom of the lake basin. In May^ "while the main water mass is still isothermal and at its winter temperature of ^'^C or less^ the edge of the lake begins to warm under the influence of a combination of factors, including warm-water run- off from the adjacent land. The nearshore boundary water warms to a temperature about S^'C (and is isothermal to the bottom) in water depth of 20 m, and extends a mile or more offshore. The edge between the boundary water and the main water mass of the lake is characterized by a sharp temperature gradient (the "thermal bar"), water color differences, marked differences in phytoplankton populations, and strong geostrophic currents along the temperature gradient. As the season progresses, the boundary water edge becomes more complex and often shows multiple bands of water color and temperature 2 to ^ miles or even more offshore, parallel- ing the perimeter of the basin. The period of 28 May to 5 June 1967^ when this process operated, has been described for Lake Michigan by Stoermer (I968) and by Noble and Anderson (1968). This same process has been well documented for Lake Ontario and Lake Huron by Rodgers (1965,1966). In the Great Lakes area, the wind system is very much modified by the con- tinuous presence of a high pressure in the spring- summer period and a low pres- sure in the fall -winter period due to the meteorological effect of this world ^s greatest freshwater masso The wind observed in the lake is quite different from that of nearby land stations o Both the magnitude and the direction of wind are frequently changing, seldom consistent over a long period of time^ The time- averaged net energy put into the lake by wind stress may not appear as a simple wind-current relationship. This is especially obvious during the lake warmup season. In the spring, when the lake water is still relatively much colder than the already warmed adjacent land, the air mass above the lake is frequently in a stable condition. Wiresonde cross- sections in Lake Michigan have demonstrated the intense springtime inversion prevalent before the warming of the water sur- face (Bellaire I965) . At the Wisconsin shore, it has been observed that water remains smooth even beneath offshore winds of greater than 10 knots (at 50 ft). However, as far as the thermal energy connected with temperature difference be- tween land and lake and air and water are concerned, temperature plays a much more constant role than the varying wind all the year around. During the lake warmup season, temperature gradients, either horizontal or vertical, show per- sistence, which, in turn, suggests that the thermal energy input may dominate the lake circulation dynamics, 59 Noble (1967) has indicated that the circulation dynamics of Lake Michigan are probably dominated by geostrophic forces rather than the wind stress during the summer season. It is the purpose of the present study to model Lake Michi- gan, considering the geostrophic thermal gradient force as the prime source for the lake current generation, and to investigate the thermal circulation pattern of the lake as a comparison with the general current reported by others (e.g., Ayers et al. 1958). The study has not reached the final stage yet. The model may be refined by some modification much closer to the real simulation of the lake. The preliminary result shows that the thermal circulation pattern in Lake Michigan is in very good agreement with general current data reported earlier (Ayers et al. 1958; Rodgers 1965) as well as our own investigations. MODELING CONSIDERATIONS AND ASSUMPTIONS In the geophysical fluid model of large scale motion, Coriolis effect must be taken into consideration. During the spring warming season, a sharp hori- zontal temperature gradient always exists within a few miles from shore around the perimeter of the lake. In a stable atmosphere, when winds are constantly shifting above the lake, the thermal gradient force comes to equilibrium with the Coriolis force. Generally speaking, only the geostrophic thermal gradient system is apparently in operation during the spring warmup season. Equilibrium between the temperature gradient force and the Coriolis acceleration produces a component of horizontal velocity perpendicular to the direction of the imposed thermal gradient in a rotating system. This is the geostrophic thermal gradient relationship which is considered dominating. In Lake Michigan, the length (about 500 km) is much greater than the width (mean about 120 km). The wide Straits of Mackinac connect the northern end of Lake Michigan to an even larger lake — Lake Huron. In the Lake Michigan model, to the lower order of approximation, the flow field, away from the north and the south ends, can be considered as zonally independent. The Coriolis force varies with geographical latitude. In Lake Michigan, its magnitude ranges from 0.^88 x 10 rad/sec in the southern end to 0.525 x 10 rad/sec in the northern end. The variation from south end to north end is small. The Beta- effect in the Lake Michigan model can be neglected. An east-west cross section of the central part of Lake Michigan is shown in Figure 5-1. The deepest point in Lake Michigan is 28I m in the northern half of the lake. The southern portion has a maximum depth of about I50 m. In the Lake Michigan model, with a width of 120 km and depth of I50 m, a flat- bottomed trapezoidal cross-section as shown in Figure 5-2 is not too far from reality. Since the dimensional reference length, the width of the lake, is small compared to the radius of the earth, the centrifugal acceleration can be neglected. 60 o •H -P u o ft u -p 0) o -p Ch o i:3 O •H -P O <D W I m o o -p m a) is +3 , ra n 01 01 W w •r) ^ • t) H •H 1 ^ ir\ a) (1) ,y ^H 01 g, h-i •H Ch tq o 61 T D 1 EAST Figure 5-2. Cross-section of a long, symmetrical trapezoidal Lake Michigan model. During the spring warming season, though the temperature structure in Lake Michigan is in a continual changing semi- steady state, the rate of heating in- creases with only some slight change in magnitude. Therefore the thermal body force is a seasonal function varying slowly with time. In modelling the lake with respect to the warmup season, all variables can be considered as time independent without losing the good approximation. The frictional forces are important near all boundaries. In the flow field of a turbulent nature as in Lake Michigan, the frictional forces are expressed as the space rate of change of the products of eddy viscosity and Reynolds stresses in analogy to the laminar flow. The average horizontal eddy viscosity in Lake Michigan is the order of 10^-10^ cm^/sec and the vertical eddy vis- cosity (1-10^ cm^/sec) is much smaller. Eddy diffusivity may be of the same order as eddy viscosity but it is a bit smaller in general. In our preliminary model 10^ cm^/sec is used as the average value for eddy viscosity and eddy diffusivity in both horizontal and vertical direction for ease in management. In the future refinement, separate values for horizontal and vertical viscosity are necessary in the Lake Michigan model. Eddy viscosity and eddy diffusivity of the fluid depend not only on temperature (such as the molecular viscosity), but change from place to place and time to time, depending on the local charac- ter and dimension of the flow. The overall average values of the order of mag- nitude of eddy viscosity and eddy diffusivity can be treated as "constant" for application to the geophysical fluid problem (see section III of this report). In general, though a weak thermocline exists during the later portion of the warming-up season, the lake is homogeneous in this particular period. In our preliminary model, we shall treat the fluid as homogeneous. Boussinesq approximations are used in the formulation, and the equation of state is a linear relation between the density and temperature as p' - P^[l - a(T' - T^)], 62 where a. and T are specified in Table 5-1. Furthermore, the dissipation func- tion is neglected in the heat equation. TABLE 5-l» Nomenclatures and nondimensional parameters. Symbol Definition X Y Z q. p V AAA Ot A Aq a L D AT Horizontal east-west coordinate, meridional direction Symmetrical north- south coordinate, zonal direction Vertical coordinate, negative is in the direction of gravity Velocity with components (u^v,w) Deviation of the pressure from hydrostatic pressure Temperature Constant lowest temperature near the bottom at the center of the lake Two-dimensional gradient operator, i -r- + k -r- ox oz (N.B. The above symbols primed have dimensions, unprimed are dimen- Two-dimensional Laplace operator, (N.B. The sionless. ) Unit vector in the (X,Y,Z) directions Acceleration due to gravity Coefficient of thermal expansion; always positive except below ^°C Kinematic eddy viscosity Thermometric eddy diffusivity Angular velocity of the vertical component of earth ^s rotation in Lake Michigan Width of the lake^ reference horizontal length dimension Depth of the lake, vertical dimension Total horizontal temperature difference between the edge and the central portion 65 TABLE 5-1 (Continued) Symbol Definition ^ Stream function for the velocities P Fluid density at temperature T 7 Depth to width ratio, 7 = D/L also nondimensional surface depth \ Tangent of the slope angle with the vertical X = tan (see Figure 5-5) C Dimension of velocity C = (XgDAT/2QL p Nondimensional thermal Rossby number, p = agDAT/(2n)^ L^ G Nondimensional Ekman number, e = A/2nL^ Eddy Prandtl number, a = A/Aq MATHEMATICAL FORMULATION With all the considerations in the previous paragraphs, let us assume that the fluid is constrained in a long symmetrical, trapezoidal, cross- section model of the lake. It is also subjected to two oppositely imposed horizontal temperature gradients toward the center from both sides of the lake. The lake center is taken as the center of the Cartesian coordinates (Figure 5-2). The lake is rotating with angular velocity Q, (component of the earth's rotation nor- mal to the lake surface) about its vertical axis, which is anti-parallel to the gravitational force. The horizontal surface of the lake is considered free with no wind stresses. The mean flow is assumed to be zonal, thus the basic geostrophic thermal gradient relation is satisfied. All specifications of nomen- clatures are summarized in Table. 5-l» The governing equations of the lake model are then as follows: q^.V^q* + 2nk X q* + i V^p» - ag(T^ - T )k = AV»^J% (l) AqV^^T' - q*-V'T^ = 0, (2) V»-q^ = 0, (5) 6h which are respectively the Navier-Stokes equations of motion, the energy equa- tion, and the continuity equation of the flow field. Note that the V and V*^ are two-dimensional gradient and Laplace operators. The above equations are, by axisymmetry, independent of Y. The boundary conditions are T' = T , q' = at z' = 0, (k) T' = -(1 - cos 2jtX')AT + T , ~^ = o at z = 7 (5) 2 o oz ||i- = , q' - at x= +1 L, (6) which imply that a horizontal temperature profile of the form of an inverse cosine curve is impressed upon the surface of the fluid of the lake, and the lake is insulated at both edges with only a slight conduction between the lower portion of the lake sides to the deeper layer of the earth. These differential equations together with all boundary conditions are coupled high order and nonlinear. Analytical direct mathematical solutions are very difficult. There- fore, some approximate approaches capable of making the equations more tractible are preferred. Nondimensional Equations It is convenient in the approximate method to introduce nondimensional variables and parameters. The following unprimed nondimensional variables are defined: rri t _ rp L ^ D ^ AT ^ ^ C C^p ^ o where C has the dimension of velocity. Since in our basic assumption the mean flow field is in geostrophic thermal gradient equilibrium, then the geostrophic thermal gradient relation: 20 — = ag — dz "=" dx leads to agD c = ^^-^'Sx _ ggPAT 20, 2nL ' where D is the depth of the fluid. 65 Substitute all nondimensional variables into equations (l), (2), and (5) with the change of vertical dimension D to the same dimension as the horizontal^ L, but using 2=7, the depth to length ratio^ to denote the free surface, the governing equations become: p(q-V)q + k X q - kT + pVp = eS^q, (T) gV^T - crpq-VT = 0, (8) V-q = 0, (9) where p _ ^ ^ is the nondimensional thermal Rossby number as the (h) ^ ^ of Hide^s (ISSk) parameter e = p is the nondimensional Ekman number a = A/Aq is the eddy Prandtl number Both p and e are very small number and a has the order of unity. The fluid is now bounded in --<x<7: ^0<z<7 where 7=7 , and boundary conditions are q - , T = at z = 0, (lO) ■^-0 ,T=-(1- cos 2tix) at z = 7 (H) oz 2 rVr 1 q = on the side boundary and ■T: = Oatx = +— . (12) Since p is small, the ordinary perturbation method can be used. We may intro- duce a stream function by letting 66 t" = - JX Vt(x-Z) + €"M y (13) where the horizontal zonal velocity v is much greater than the meridional velocity, u, and the vertical velocity, w, which will appear in the following sections. Expand all the variables as powers of p as vn n (n) ^ n=0 "^ ^ n=0 ^ n=0 ^ where the superscript n denotes the nth order of p- expansion. Equate terms in the powers of p. The zeroth p- order equations are: W°^ +v^°^ -T^°^ =0 (lU) Z X eV-;°) - A°) =0, (15) W°^ = 0. (16) The first p- order equations are ^4(1) (1) (1) r (O) (O) (O) (O) (O) (O) (O) (O), Z X Z XZZ Z XXX X zzz X xxz (17) eW^) -^^'^ - [/°)v(°^ -^(°V(°)] = 0, (18) Z Z X Z Z eW") - CT[^^°^T^°) - ^^^°^t(°)] = 0, (IV) Z X X z and similarly for the other high p- order terms. 67 Boundary conditions are that all p- orders of velocities, and in turn, all p-orders of stream functions, should separately satisfy equations (lO), (ll), and (12); and that all p-orders of temperature should satisfy (lO) and (12) with (ll) for the zeroth p- order and T^^ = 0, (20) for other higher ith p-order temperatures, i = 1,2,3, etc. SOLUTIONS OF THE ZEROTH p- ORDER APPROXIMATION Temperature The zeroth p-order temperature can be obtained from equation (l6), which is simply the Laplace equation. The solution of equation (l6) satisfying boundary conditions (lO), (ll), and (12) is (o) Irz sinh 2nz ^ - /^ v r ^ = -L- - . ^ ^ cos 2TtxJ . (21) 2 7 smh 2tC7 The zeroth p-order temperature is purely conduction. Boundary Layer Scaling Equations (ih) and (15) are coupled and show boundary-layer characteristics because e is a very small number. In boundary- layer analysis the boundary- layer thickness should first be determined as a function of the small parameter, e. Equations (1^) and (15) imply that eW°^ + tl°^ = . (22) On the free surface and the rigid bottom we may scale the normal coordinates near the boundaries as cp :. €"^(z - 9) (25) where = on the bottom and = 7 on the surface. 68 Equation (25) substitutes into equation (22) with no stretching of the horizontal coordinate. The small parameter €, multiplying the highly differ- entiated term is important only near the boundary where the frictional force is the same order of magnitude as the other terms. This implies a = l/2. On the inclined side boundaries ( see Carrier 1955) where X is the slope of the slanted boundaries, and X = tan as shown in Figure 5-2; let I = ^""^(x + 0), (25) and where = +[2 - X(7 - z)] on each of the two side boundaries. By substituting equation (25) into equation (22), to the lowest order of (e), we obtain b = l/2. Therefore, all the boundaries have the boundary layer thickness of the order Following the procedure of boundary- layer additive type analysis, let (o) (o) (o) (o) (o) (o) ... r ^ = -^r ^ + ir ^ + ^r ^ + sr ^ + ^r . (26) where the left subscript denotes the region. The jt is dominating in the interior region and -t^^^, i = 1,2,5^^ ^^^ significant only on their respective boundaries. The various boundary- layer regions and their widths are shown in Figure 5-5' Coordinates for the regions of the lake cross- section are summarized in Table 5-2. Interior Solutions In the interior region, boundary effects are small. Equations (ik) and (15), with the first terms neglected, become: v(o) Jo) ^ ,(0) ^ ^ ^ ^ = T^ ^ and ^\|;^ ^ = I X I z 69 JJL (I) 12Y (I) indicates the Interior region (1) Indicates the surface region (2) Indicates the bottom region (3) Indicates the "eastern" boundary region (4) Indicates the "western" boundary region Figure 5-5. Boundary regimes of the Lake Michigan model, TABLE 5-2. Local coordinates in different regimes. Region Coordinat es I X ^ 1 3£. CPl = € (z • ■ 7) 2 X -1/2 cps = e z 5 ll = e .1/2 [X - 1 2 + Hy - z)] \=^ k I2 = e .1/2 [X + 1 2 - Hr - z)] Tl^-Z which yield: and v(o) l.cosh 27tz . >. n «/ \ E 2 smh 27t7 -^ ^ ^ ^ V°^ =B(x), (27) (28) where A(x), B(x) are integration constants and possibly functions of x. A(x) and B(x) can only be determined by matching up the interior solutions to the solutions obtained from the boundary- layers. 70 Boundary- Layer Solutions A. Region 1 — the upper boundary (free surface) In the upper boundary region, the stream function and the zonal velocity are, f - -j-f + i¥ ^ V - V + iv and the vertical coordinate is stretched as shown in Table 5-2. Equation (22), using the local coordinates, yields Since it^ ^ vanishes away from the upper boundary, no nonzero forcing term can be satisfied after integration of equation (29). Therefore equation (29) can be simplified as: Together with the appropriate boundary conditions, as 5.^°^ ,(o) ^i^°^ , , or equivalently, equation (30) can be solved and we obtain: (o) cpi/f2 . , I- , , i\i/^ '^^ = -GTce^-^'^ sin 2Ttx cos cpi/v2 , (32) and ( o) 1/2 Tte . ^ / cpi cpix iv^ = -eV^ ^^ sin 2^.x(cos ^ + sin ^) 71 (35) B. Region 2— the lower boundary (rigid bottom) Stretching the vertical coordinate^ equation (22) leads to an equation similar to (50). With the lower boundary conditions (o) (o) S^l;^''^ "^ ^""t = at cps = , (5U) we obtain 2t' =-€:te ^^'^ sm 2Trx(cos ^ + sin ^) , (55) and M°^ = cV^vfl- «e-^^/^ sin 2«x cos ^ . (56) Matching the boundary solutions with that of the interior, the boundary conditions in equations (5I) and (^k) also yield the two integration constants in equations (27) and (28) as A(x) = sin 2Ttx(€V2 fz - - — ,^ ^ ) , ^ ^ ' 2 sinh 2Tr7^ ^ and B(x) = -€tt sin 2Ttx . Then the interior solutions are: (o) j^ = Gjr sin 2:n:x , (57) (o) , 1/2 nr cosh 2jtz - 1 , . ^ ... _v' = (g / \r2 +•- — . , ^ ) sin 2nx . (58) I ^2 smh 21(7 C. Regions 5 and h — eastern and western rigid side boundaries Using the local stretched coordinates normal to the side boundary, equation (22), to the lower order of €, leads to 72 Oil in region 5» With the boundary conditions, So) ( o) St , , . v^ ^ = ^-— = at li = (UO) 051 we obtain. ,(°) e(l + A.^)k kii/y/i" . ^ ,1 ,, ... cosh 2Ttnj - 1 r^n/ Mi • Hi\ (o) k|i/>/2' k|i . ^ a , . ,,cosh2itTii- 1 r?r 1/2-, where k = \ /(I + X ) . By symmetry, similar solutions are obtained in region k as 4' (o) e(l + \^)k -)^iz/{2 . ^ a ,, , , , cosh2rtT]2-l r~,. kfo . k|p. ,(o)_ -k|2/V2"_„ k 4V • = -e --' ■ cos (111.) All the boundary- layer solutions are significant only on their respective local boundaries and go to zero away from their own boundaries. From equation (26), the general solution of the zeroth p- order approximation is: 75 .^^°^ = e[n sin 2^[-l . e^'^/^^ cos ^^ . e"^/'^ (cos -;#= . sin ^) ] _-.[..l/2-X(.-.)„ ^ (^^, J°) • o • r cosh 2itz - 1 ^ 1/2 ,r- rt -J^e^ '^^ rrr -z/\f2e z , ^ = ^^^ "'^f 2 sinh 2.r """^ ^^W^^ -^^^ °°^-;|??] ^ {sin 2n[i - A.(7 - ,)](<^°^^2^^- 1 + ^1/2^), ^^-k[x4-l/2-X(r-z)]/v[2r 2 2 sinh 2tC7 - cos kfx + 1/2 - X(y - z) 1 k[x-l/2+\(r-z) ]/\[2r kfx - l/2 + A.(y - z) 1 cos ^ - e cos ^ ) ' and equation (21) is the zeroth p-order temperature. RESULTS OF THE ZEROTH p- ORDER SOLUTION For a small thermal Rossby number^ p^ and a high Taylor number^ equivalent to € ^ our solutions are within the range of lower symmetrical regime in a horizontal differentially heated rotating fluid (Fowlis and Hide I965). The governing equations are linearized by ordinary, perturbation expansion with p. Solutions are obtained by further singular perturbation expansion for each order of p in performing boundary- layer analysis with respect to the Ekman number €. Temperature The zeroth p-order temperature in equation (21) is obtained from heat con- duction, neglecting the convectional and the frictional effects. The temperature distribution is a linear function of depth peaked on the surface edge with some horizontal variation minimized at the central portion of the lake. Figure 5-^4- presents the horizontal temperature distributions of available data at the Great Lakes Research Division from 1965 to 1968 during the warm-up season. Figure 5-5 shows the vertical temperature distributions during the same period of time in the central portion of Lake Michigan. Generally speaking, both the horizontal and the vertical temperature distributions show the same pattern as the solution indicated. The solution is in good agreement with the mean synoptic temperature measurement, especially in the interior region. On the boundaries, either hori- zontal or vertical, the agreement is not as good as the interior. On the upper ::::!:-:9 12.0 10.0 10.0 -43** 30' JUNE 1965 GRAND HAVEN 43«00* 19.0 17.0 - 43*^30' JULY GRAND HAVEN 43° 00- ■ 43»30' AUGUST \v \ 1/'^^^^^^^ HAVEN ' ' \ V yy^y^y 43*»oo' (a) Figure '^-k. Surface temperature patterns during spring warming of Lake Michigan from available data, I95U-68. 75 87O00' 4.0 43*»30' MAY 1966 GRAND HAVEN ::•::.' 43« oo' 8.0 .6.0 V. ^^^ _ — — / 7.0 \ / / / / 1 \ 7.0 \ / \ 43** 30' JUNE GRAND HAVEN . 43«00' 43'* 30' AUGUST GRAND HAVEN X 43** 00' 20.0 (b) Figure 5-^ (Continued), 16 ■ 43° 30' MAY 1967 GRAND HAVEN 43° 00' ■ 43° 30' mi I ' f / \ I \ \wmm JUNE il ' i 1 1 W X'>7 / . / /I5 0)^„ I }i50^^^^^^^ HAVEN Figure 5-1^ (Continued). 77 6.0 4.0 ► •"•"•'•"• « •••••« • ••••••• aV ■ x-:-:-:v:vy * 14.0 \ • •••••• • \ ' *\j,\j % 12.0 87*»00' iig^im::.. MAY 1968 43*' 30' GRAND HAVEN 43** 00' 43** 30' Figure ^-k (Concluded). 78 CO o CO LiJ UJ < a o •H -P ^ O Ph H (d ^ •P S <u C) • 00 «H VO o 1 \r\ 0) M3 ^ ON H 4:> o •\ g t -p •H QQ Jh (U ft b bO -P d cd ^ g 0) i Q) +3 bO a H •H cd O 1 •H TJ 1^ 0) § > ttO :d k o g sj a; 9UJ 79 CD Z o CD LiJ CO Z) llJ < 80 MILWAUKEE MUSKEGON MAY MUSKEGON 16° I7*»_ IP- 17° JUNE 19° 18* 1967 I9« 19° 2»-l8° AUGUST Figure 5-5 (Continued). 81 00 CD CD o UJ o 00 o g 5 ^ J^ UJ 1 ^^ 1 1 S SJ949UJ 82 horizontal boundary, the real temperature decrease with depth is more rapid than the relation predicted by the solution. On the lower bottom layer, the real temperature decrease is much less. On the two slant side boundaries, real temperature gradients exist quite close to the shore, and in the central portion of the lake the temperature change is very small. The temperature variation is much greater than the solution yields. All these observations may indicate: (1) the impressed surface temperature is not an exact inverse cosine profile, and (2) the temperature is modified by convection as a result of the coupling of the temperature with the inertial terms in the equation of motion, particularly on the boundary layers. For mathematical simplicity, by smoothing the tempera- ture gradient, the inverse cosine temperature distribution is not too far from reality. As for the importance of convection in determining the temperature, higher p-order approximations should be attempted. Since the zeroth .p-order velocities have boundary- layer characteristics, and these velocities together with the zeroth p-order temperature gradients are forcing the first p-order temperature, it is expected that the higher p-order temperature will show boundary-layer characteristics and modification of the overall temperature dis- tributions. Velocity From the basic realistic assumption of the geostrophic thermal gradient balance, neglecting variable surface wind stresses, the zonal velocity dominates the overall circulation of the whole lake as expected. In the interior region, the pressure force balances the Coriolis force due to the zonal velocity together with the thermally induced body force. This basic assumption of geostrophic thermal gradient relationship (the balance of the vertical zonal velocity gradient and the horizontal temperature gradient) must be maintained. The interior solu- tion of the zonal velocity is shown in equation (38). It has the maximum magni- tude around one-fourth the width of the lake from the shore edge near the sur- face and approaches minimum as x approaches zero and half indicating the middle portion and both edges of the lake. During this warm-up season, on the eastern portion of the lake where x is positive, there is a northward zonal current. On the western portion of the lake, where x is negative in the horizontal coor- dinate of the lake model, there is a southward current. Both the horizontal and the vertical velocities are much smaller (in the order of e or ei/s) than the zonal velocity. From the interior stream function equation (57), there is no horizontal velocity (u = SVSz) in the interior region. The vertical velocity, w = - d\|f/Sx, leads to a broad but very small vertical circulation in magnitude in the interior region given as (o) P ^w' ' = -€ 2Tr cos 2Trx ( J4.7) with a relatively strong sinking in the middle of the lake, and slight updwell- ing along the edges. 85 In region 1, that is^ in the free surface layer^ the zonal velocity is strengthened by adding to the interior zonal velocity a boundary zonal velocity of order (g-*-/^) in magnitude. The direction of the general zonal flow is the same as that of the interior, which shows a cyclonic circulation. The meridional velocity in the upper boundary is different from the interior where no east-west horizontal component develops. On the surface layer^ there is a meridional velocity ..(o) St -1/2 ^ 1/2 ^ ^i/{^ . ^ / 9i . ^i\ ^^' = -^- = € /^ -L_ = e/^-:^ e^^/^ sin 2^(cos ^ - sm ^) (^8) which flows from shore boundaries due east and due west^ toward the central portion of the lake. The vertical velocity in region 1 is much smaller than the horizontal velocity as the continuity equation implies. The vertical velo- city is: (o) (o) (o) (o) ^2 ^1/^2 cpi ^ .,^. W^ ^ = W^* + iW^ I^ ^ + G 2Tr e^^^"* cos ^ cos 27tK (1+9) which is the continuation of the interior velocity with only a slight increase in intensity. In region 2, the bottom layer, there is a meridional flow of an equal or- der of magnitude but opposite in direction to that of the surface, with a ver- tical velocity similar to that of the surface. The zonal velocity in the bottom layer is much smaller, as expected. It can be noticed that the residual zonal velocity near the bottom is the difference of the interior and the boundary zonal velocities. In the side boundaries, regions 5 sind U, the local coordinates (|^r|) are stretched normal to the inclined surfaces. From equation (^l) the meridional velocity in region 5, to the order of e, is: / -1/2 . S S V ,( o) ,(0) _ ^t^°^ that is (0) -eV^ k^(i + x2)^k|i. . k^i/>/2' . k^i . ^ , v.cosh 21^111- ^ rT\ )e sxn -f sxn 2^7 " ri3)( , ^^^^\^ - <^) (50) &k This meridional current flows away from the shore with a moderately high speed in the side boundary as |i is negative when X is positive, and vanishes in the interior region. By symmetry, a similar flow away from the west coast toward the center of the lake exists in region k. The local vertical velocity in the eastern alongshore region, to the order of e, is ,w (o) ax = ._ei/2 k^d + x^) e^^^/ ^ sin^ sin 2TtX(r Tll)( cosh 2Trr]T - 1 2 sinh 2ny (51) which represents sinking in the side bounda ry- layer . The local vertical velo- city in the alongshore boundary region is opposite to the vertical velocity in the interior. Therefore in the alongshore boundary region^ the overall vertical velocity may be small. Similar symmetrical phenomena exist in region k. The local zonal velocity in regions 5 and k, as shown in equations (^2) and {^h), are opposite in direction to the general pattern of the zonal circu- lation extended to the boundary from the interior. In the east and west along- shore boundary region^ the current direction is dependent on the difference in magnitude of the general zonal current and the local zonal current. Therefore^ the overall zonal currents in these two alongshore boundary regions may flow either southward or northward. Theoretically^ boundary currents very close to the shore should be smaller. In summing up the foregoing current analyses^ the general horizontal zonal circulation is shown in Figure 5-6. The meridional and the vertical circulations are shown in Figure 5-7a- When the overall temperature in the lake is less than i^.°C^. the coefficient of thermal expansion of water^ a, is negative. Hence the meridional and vertical circulation patterns of the lake are changed as approx- imately shown in Figure 5-Tt>* WEST SHORE A^\ I r i J I I (1) (I) ♦ (2) f f.l 1(3) EAST SHORE Figure 5-6. General horizontal zonal currents. 85 WEST -:^ * (i)* * 5= 5r-T=2 i i i(I)l I i f I t * - ^ *. *(2) * _^ _^ ^ ^ v. V> V^ v^ v^ v^ v^ *^ ^^ V^^ V a) . For a positive EAST ^ -^ *- — —-.- (1) * - WEST i t » t I'^'Im ^ ^ ^ ^ (2) . . V^ V. v^ V. b) • For a negative Figure 5-7. Meridional and vertical circulation. EAST Since we have not analyzed in detail the magnitude of each component of velocity. Figures 5-6 and 5-7 only roughly show the possible general circula- tion pattern. A synoptic surface current computed from dynamic heights in June 1955, reproduced from Ayers et al. (1958), is shown in Figure 5-8. At the central portion of the lake, the general circulation pattern in Figure 5-6 is in good agreement with that of Figure 5-8. A series of drogue studies for surface current measurement are shown in Figure 5-9 for the eastern shore area, and in Figure 5-10 for the western shore area. Tables 5-5 ^^^ 5-^ give the elapsed time for the drogue runs and the average drogue speeds for the east- shore and west-shore experiments, respectively. Generally, drogues drift toward the north near the eastern shore and toward the south near the western shore, which is in good agreement with the theoretical solutions. Winds are also in- dicated in Figures 5-9 and 5-10. In some cases, currents were flowing directly against the winds. In studying the thermal bar phenomena (the temperature gra- dient zone) in the Great Lakes, Rodgers (1965) indicated that the circulation associated with the thermal bar should be in the pattern as shown in Figure 5-11 (Figure 7 of Rodgers* paper). When compared to our theoretical current pattern in a vertical section shown in Figure 5-7ti the major characters of these features agree well in general. Of course, our model does imply that this kind of circu- lation pattern exists only when the temperature in the thermal bar is lower than 4°C. It is expected that when the sharp temperature gradient exists around the k°C zone, this kind of circulation pattern is more augmented by the acceleration of the vertical velocity in the boundary by the negatively buoyant force. 86 Figure 5-8. Surface currents in Lake Michigan, 29 June 1955 (after Ayers et al. 1958)- 87 86°30' 43°10' (a) / \ B. 86^20 (b) 86°30' 43**00' / K. L. M. V 86°20» 43^00* Figure 5-9 • Drogue studies for surface current mea,surements near eastern shore, (a) 2 May I968 (wind N, 10 kts?). (b) 2k June 1968 (wind SW 15 kts). 88 (c) 66°30» 43°10» Date TABLE 5-3 Drogue Studies on the East Side. Wind Drogue Station Start Finish cm/sec. 2 May '68 N-10 24 June SW-15 A. B. C. D. E. F. G. H. I. J. K. L. 0950 1003 1009 1015 1020 1026 1031 1032 1027 1022 1017 1012 1504 1514 1522 1534 1545 1554 1559 1502 1520 1548 1602 1612 7.14 5.80 5.80 3.12 6.24 6.24 4.01 8.94 5.80 7.70 10.26 8.96 30 July SE-20 M. 1007 1623 9.41 N. 1002 1637 9.41 0. 1005 1328 28.9 P. 0952 1300 19.5 Q. Lost R. 0931 1220 10.1 S. 0918 1235 12.7 (c) 50 July 1968 (wind SxE, 20 kts?) Figure 5-9 (Concluded). 89 43*'00 Milwouke* (a) CALM E. F. T^ ^ TABLE 5-4 Drogue Stations on the West Side. Wind Drogue station Start Finish cm/sec . 10 May '68 Calm Miles 27 June 87^50' NNW-6 43*^00' \ A. 0800 1506 3.12 B. 0815 1520 4.73 C. 0821 15 32 8.02 D. 0837 1548 4.91 E. 0848 1604 7.58 F. 0859 1620 5.36 G. 0912 1635 4.01 H. 1245 1848 4.91 I. 1251 1820 9.86 J. 1301 1813 10.74 K. 1308 1805 7.14 L. 1314 1856 4.91 M. 1321 1840 8.47 N. 1327 1831 7.14 (b) 43° Figure 5-10. Drogue studies for siirface current measurements near western shore, (a) 10 May I968 (wind slight), (b) 27 June I968 (wind NxW, 6 kts). 90 WARM^N^ Figure 5-11. Circulation associated with "Thermal Bar" (after Rodgers I965). THE FIRST p- ORDER SOLUTIONS As we have already mentioned, the first p- order governing equations of the model are equations (1?)^ (18)^ and (19). The forcing terms which are the products of the zeroth p-order temperature gradients and the velocities (as shown in equation (19))^ are of the order unity in the interior region, of the order (e"""/^) for the individual local boundary- layer regions at their re- spective boundaries, and with another forcing term of the order unity due to the overall interior influence. Since the forcing terms are of the boundary- layer type, it is expected that the first p- order temperature will have strong boundary- layer characteristics. In each of the boundary regions, the first p- order temperature can be separated into two parts; one part is defined by the overall existing zeroth order interior stream function (^t^^^), and the other part is due to the local stream function which is important only in its own boundary region. In this way the boundary- layer portion which is forced by the terms of order of (e ^/^) can be treated with boundary- layer analysis. After applying the respective local boundary conditions, solutions are obtained for all the first p- order temperatures influenced by local zeroth p- order stream functions. Then at all the new boundary conditions, solutions are obtained by applying the first p-order solution which has the boundary- layer characteristic as a perturbation of that part of the temperature under the Laplace operator (which has no boundary- layer characteristic). In the interior region the tem- perature is separated as a general solution and a particular solution, with the forcing terms of the order of unity operating only on the particular solution. After transferring the nonhomogeneous terms as a particular solution into the new boundary condition and matching with the boundary solutions to the order of unity, the overall first p-order temperature can be solved. The result of the first p-order temperature turns out to be of the order of unity and there are some modified variations of the order ( e^/^) in the boundary regions. 91 After obtaining the first p-order temperature, the first p-order stream function and the zonal velocity can be obtained from the coupled equations (I7) and (18) following the method used in the zeroth p-order solutions, and with more complexity. CONCLUSIONS After the winter mixing, water in Lake Michigan is almost homogeneous. The spring warming season begins during late April or early May, ' Beginning at this time the lake starts to warm from the shallow water along the shore progressively outward to the main body of water which is still at k^'C or lower. During this period, the land and air are much warmer than the water. A stable situation usually exists in the atmosphere above the lake. It is expected that much less energy is transferred from the wind to the lake water. A sharp tem- perature gradient always exists within a few miles along the shore which leads to a balanced geostrophic thermal gradient model of Lake Michigan in a simpli- fied long trapezoidal shape. A double perturbation method is used in obtaining approximate solutions from these coupled nonlinear momentum and energy equations by expanding first with the small thermal Rossby number, p, and then, to each order of p- using singular perturbation, with respect to the small Ekman number, €. All boundaries are found to have the same order of boundary- layer thickness in the order of ( € ^f). The zeroth p-order temperature is found to be due en- tirely to conduction. The first p-order temperature shows that convection is important and it has boundary- layer variations. All the velocities do have boundary- layer characteristics. The general zeroth p-order horizontal zonal velocity shows a cyclonic circulation pattern with variations in the along- shore boundaries. When the overall temperature of the lake is still lower than k°C, the horizontal meridional velocities and the vertical velocities show a broad upwelling of a small order of magnitude in the middle portion of the lake and a sinking zone along the edge of alongshore boundary- layers coupled with the divergent horizontal surface meridional current and the convergent bottom cur- rent. This forms two gigantic but weak counter-rotating cellular motions. Two horizontal meridional currents flowing outward from their respective shores in the two alongshore boundary regions to joint the thermal bar convergences sug- gest the possible existence of another upwelling very close to the shore bot- toms. Though detailed characteristics have not been wholly demonstrated, the general circulation pattern compared with results obtained by Ayers et al. (1958) sind Rodgers (1965) indicate very good agreement. From a physical point of view, in the spring warming season the shallow water alongshore region re- ceives a great amount of heat energy from the air and the land. The temperature in these regions increases very rapidly while that of the main body of water does not respond as quickly and remains at the temperature of maximum density. This temperature gradient induces an inclined surface water slope tilted up toward the land, which induces a "Hadley cell" (1755)- The thermal current is originally flowing outward from shore but is deflected. by the Coriolis force which results the northward or southward zonal current. Due to the relatively small dimensions, the boundary effects are more important in lakes than in 92 oceans. Though we have not analyzed the numerical magnitude of all variables^ the general circulation pattern has demonstrated all the main expected results, FUTURE PIAN Our preliminary model is too rough to obtain detailed characteristics of the lake circulation either in magnitude or in flow pattern. It is my hope to be able to further refine the model both in dimension and magnitude of some var- iables in similation of the prototype of the lake. For example^ the horizontal dimension and the vertical dimension should be scaled separately; the vertical velocity should have a different velocity scale than the horizontal one; the horizontal eddy viscosity and the vertical eddy viscosity should have different mean values. As in perturbation procedures, higher orders of approximation are very important. The detailed analyses of the first g-order solutions and the estimation of approximate error from the second P-order equations are scheduled to be carried out. The combined final solutions should be programmed into com- puter calculations and the results should be compared with all available field data including some of the current meter buoy stations of the Great Lakes- Illinois River Basins project in 196^. Furthermore, I also intend to consider the thermal body force as a slowly time dependent function and to extend the application of the model for the entire year. 95 References AYERS, J* C- D. E. CHAHDLER, G. H. LAUFP^ C. F. POWERS, and E. B. HENSEN. 1958. Currents and water masses of Lake Michigan. University of Michigan, Great Lakes Res. Inst. Pub. No« 5, I69 p. BELLAIRE^ F. R. I965. The modification of warm air moving over cold water. Proc. 8th Conf., Great Lakes Res., University of Michigan, Great Lakes Res. Div. Pub. No. 15, pp. 2^19-256. CARRIER, G. F. 1955* Boundary layer problems in applied mechanics. Adv. in Appl. Mech., 3:1-19. FOWLIS, W. W. and R. HIDE. 1965. Thermal convection in a rotating annulus of liquid: effect of viscosity on the transition between axi symmetric and non- axisymmetric flow regime. J. Atmos. Sci. 22:5Ul. HADLEY, G. 1755. Concerning the cause of the general trade-winds. Philo Trans. Roy. Soc.j, London, 59:58-65. HIDE, R. 196U. Experiments on baroclinic instability and "vacillation": a review. Trans. Am^ Geophys. Union ^5: 1^+2-1^+5* NOBLE, V. E. and R. ANDERSON. I968. Temperature and current structure in the Grand Haven, Michigan, vicinity during thermal bar conditions. In press, Proc. 11th Conf. on Great Lakes Res.^, Internat. Assoc, for Great Lakes Res. NOBLE, V. E. 1967. Evidences of geostrophically defined circulation in Lake Michigan. Proc. 10th Conf. Great Lakes Res., Internat. Assoc, for Great Lakes Res., po 289-298. RODGERS, G. R. 1965. The thermal bar on the laurentian Great Lakes. Proc. 8th Conf. on Great Lakes Res., University of Michigan Great Lakes Res, Div. Pub. No. 15:558-565. _____^ . 1966. The thermal bar in Lake Ontario, spring 1965 and winter 1965-66. Proc. 9th Conf. on Great Lakes Res., University of Michigan, Great Lakes Res. Div. Pub. No. 15:569-57^- STOERMER, E. F. I968. Near shore phytoplankton populations in the Grand Haven, Michigan, vicinity during thermal bar conditions. Proc. 11th Conf. Great Lakes Res., Internat. Assoc. Great Lakes Res. In press. 9h